## PRELIMINARIES 
# Need to run this code in R 64
# Need to install WinBUGS (http://www.mrc-bsu.cam.ac.uk/software/bugs/, and follow the instructions)
# and choose the appropriate path for the bugs.directory. Note that, for Windows 10, some problems
# have been reported when installing WinBUGS in "Program Files"; it is better to install WinBUGS in 
# an alternative folder or in the Desktop. 
# Also need to install the following R packages: ggplot2, MCMCpack, mvtnorm, spam, R2WinBUGS

options(warn=-1)
setwd("//isad.isadroot.ex.ac.uk/UOE/User/Desktop/PSRM - Katz and Levin (2016)/Replication")

library(ggplot2)
library(MCMCpack)
library(mvtnorm)
library(spam)
library(R2WinBUGS)

rnd = function(x, n) {
  posneg = sign(x)
  z = abs(x)*10^n
  z = z + 0.5
  z = trunc(z)
  z = z/10^n
  z*posneg
}


################################## FIGURE 1 #####################################################################################################################################################################################################
df<-read.csv("district.1986.2014.csv")
unique.state<-unique(df$State)

TREs<-c()

for (i in 1:length(unique.state)) {
TREs<-rbind(TREs, c(as.character(unique.state[i]),
quantile(df$TRE[df$State==unique.state[i]], 0.25),
quantile(df$TRE[df$State==unique.state[i]], 0.75),
median(df$TRE[df$State==unique.state[i]]),
quantile(df$TRE[df$State==unique.state[i]], 0.25),
quantile(df$TRE[df$State==unique.state[i]], 0.75)
)
)
}

TREs<-data.frame(TREs)
names(TREs)<-c("State", "lower", "upper", "middle", "ymin", "ymax")

for (k in 2:ncol(TREs)) {
TREs[,k]<-as.numeric(levels(TREs[,k]))[TREs[,k]]
}

p.1 <- ggplot(TREs, aes(x=State, lower=lower, upper=upper, middle=middle, ymin=ymin, ymax=ymax)) + 
geom_boxplot(stat="identity", outlier.shape=NA, fill="gray80")+scale_y_continuous(limits=c(0,3.5))+
theme_bw()+ylab("TREs' clearance rates")
p.1<-p.1+theme(axis.title.y=element_text(vjust=1.25),
axis.title.x=element_text(vjust=-.35))+geom_hline(aes(yintercept=mean(df$TRE)),linetype="dashed")

# Figure 1 - display #
    p.1  

# Save figure 1
  jpeg("Figure1.jpg")
   p.1
   dev.off()

## Change in TREs over time
(mean(by(df$TRE, df$Year, mean)[6:8])/
mean(by(df$TRE, df$Year, mean)[1:3]))-1


################################ FIGURE 2 #####################################################################################################################################################################################################
df.depvars<-subset(df, select=c("Invalid", "Absenteeism", "Year"))
df.depvars<-df.depvars[order(df.depvars$Year),]
df.invalid<-subset(df.depvars, select=c("Year", "Invalid"))
df.invalid<-cbind(df.invalid, "Invalid Votes")
names(df.invalid)[2:3]<-c("Value","Votes")

df.absenteeism<-subset(df.depvars, select=c("Year", "Absenteeism"))
df.absenteeism<-cbind(df.absenteeism, "Electoral Absenteeism")
names(df.absenteeism)[2:3]<-c("Value", "Votes")

df.abstention<-rbind(df.absenteeism, df.invalid)
df.abstention[,2]<-df.abstention[,2]*100


p.2<-ggplot(df.abstention, aes(x =Year, y = Value)) + 
            geom_point() + 
            facet_wrap( ~ Votes, scales="free") + 
            geom_smooth(colour="gray40", size=1.5, color="gray50", 
		fill="gray20", level=0.9) + 
		scale_x_continuous(breaks=c(1986, 1994, 2004, 2014))+
	xlab("Election")+ylab("%")+
	theme_bw()+theme(strip.text.x = element_text(size = 12, face="bold"))


# Figure 2 - display #
   p.2

# Save figure 2
  jpeg("Figure2.jpg")
   p.2
   dev.off()

################################# TABLE 1 ######################################################################################################################################################################################################
################### MCMC algorithm for individual-level model ###############################################################################################################################################################################

# Need to run MCMC algorithm first

df.individual<-read.csv("individual.2002.2010.csv")

state<-as.numeric(factor(df.individual$State), levels=unique(df.individual$State))
election<-as.numeric(factor(df.individual$Year), levels=unique(df.individual$Year))
num.districts<-length(unique(state))
num.elections<-length(unique(election))

Y<-cbind(df.individual$Valid, df.individual$Invalid,df.individual$Absenteeism)
y<-c(sapply(1:nrow(Y), function(i) which(Y[i,]==1)))
n<-length(y)          # Number of obs
K<-3			    # Number of outcomes

Matrix.Election<-as.spam(matrix(0,n,num.elections))           
for (i in 1:num.elections) Matrix.Election[election==i,i]<-1
Matrix.District<-as.spam(matrix(0,n,num.districts))
for (i in 1:num.districts) Matrix.District[state==i,i]<-1

X<-cbind(1,   df.individual$Income, df.individual$University,
		df.individual$Secondary,
		df.individual$Political.knowledge, df.individual$Urban,
		(df.individual$Candidates-mean(df.individual$Candidates))/(2*sd(df.individual$Candidates)), 
		df.individual$Political.inefficacy, 
		 df.individual$No.representation, 
		df.individual$Dissatisfaction.democracy, 
		(df.individual$Growth-mean(df.individual$Growth))/(2*sd(df.individual$Growth)),	
		(log(df.individual$Inflation)-mean(log(df.individual$Inflation)))/(2*sd(log(df.individual$Inflation))),
		df.individual$TRE, 
		df.individual$Illiterates, 
		df.individual$Young,
		df.individual$Seniors, 
		(df.individual$Competitiveness-mean(df.individual$Competitiveness))/(2*sd(df.individual$Competitiveness)),
		(df.individual$Age-mean(df.individual$Age))/(2*sd(df.individual$Age)),
		 df.individual$Female, df.individual$Married,
		df.individual$Religion, df.individual$Union, 
		df.individual$PT, df.individual$PSDB, df.individual$PMDB
		)
	  
colnames(X)<-c("Intercept", "Income", "University", "Secondary",
		"Pol.Knowledge",
			"Urban", "Candidates", "Political Inefficacy", "No Political Representation",
		"Dissatisfaction w. Democracy", "Growth", "Inflation",
		"Clearance Rate", "Illiterates", "Young", "Seniors", 
		"Competitiveness", 
	"Age", "Female", "Married", "Religion", "Union",
	 "PT", "PSDB", "PMDB")

# Iterations &  burn-in
 nsim<-2000000				# Number of MCMC Iterations
 thin<-150                          # Thinning Interval (should be a divisor of nsim/2)
 burn<-500000                       # Burn-in
 lastit<-(nsim-burn)/thin	      # Last stored value

# Priors
set.seed(112215)
mu0<-rep(0,ncol(X))
Tau0<-diag(1,ncol(X))
Sigma0<-solve(Tau0)

# Store Values
 Beta2<-Beta3<-array(0,c(lastit,ncol(X),3)) 
 Sigma.election<-array(0, c(lastit,4,3))
 Sigma.district<-array(0, c(lastit,4,3))
 Reffect.district<-array(0,c(lastit,num.districts*2,3))
 Reffect.election<-array(0,c(lastit,num.elections*2,3))

colnames(Beta2)<-colnames(X)
colnames(Beta3)<-colnames(X)

# Start of MCMC algorithm
for (ch in 1:3) {

# Inits
  b2<-rep(0,ncol(X))
  b3<-rep(0,ncol(X))
  cov2<-cov3<-0.01*diag(ncol(X))
  reffect.district<-matrix(0, num.districts,2) 
  sigma.district<-diag(2)
  reffect.election<-matrix(0, num.elections,2) 
  sigma.election<-diag(2)

 A<-0

for (i in 1:nsim) {
 # Calculate Likelihood based on old value of b2 and b3
  eta2<-X%*%b2+reffect.district[state,1]+reffect.election[election,1]
  eta3<-X%*%b3+reffect.district[state,2]+reffect.election[election,2]

  p1<-1/(1+exp(eta2)+exp(eta3))
  p2<-p1*exp(eta2)
  p3<-p1*exp(eta3)
  lold<-sum(log(p1)*(y==1)+log(p2)*(y==2)+log(p3)*(y==3))

 # Draw Candidates
  b2new<-b2 + rmvt(1,sigma=.01*cov2,3)  # Tuning parameter to reach appropriate acceptance rate	
  b3new<-b3 + rmvt(1,sigma=.01*cov3,3)   	

 eta2<-X%*%c(b2new)+reffect.district[state,1]+reffect.election[election,1]
  eta3<-X%*%c(b3new)+reffect.district[state,2]+reffect.election[election,2]
  p1<-1/(1+exp(eta2)+exp(eta3))
  p2<-p1*exp(eta2)
  p3<-p1*exp(eta3)
  lnew<-sum(log(p1)*(y==1)+log(p2)*(y==2)+log(p3)*(y==3))


# Acceptance prob on log scale
   r<-lnew+dmvnorm(b2new,mu0,Sigma0,log=T)+dmvnorm(b3new,mu0,Sigma0,log=T)-	
	(lold+dmvnorm(b2,mu0,Sigma0,log=T)+dmvnorm(b3,mu0,Sigma0,log=T))
   if(log(runif(1))<r){
    b2<-c(b2new)
    b3<-c(b3new)  
    A<-A+1
    }

###################################
# Update District reffects using MH #
 ##################################
 reffect.district.new<-reffect.district+rmvt(num.districts,sigma=0.001*diag(2),3) 
  eta2.old<-X%*%b2+reffect.district[state,1]+reffect.election[election,1]
  eta3.old<-X%*%b3+reffect.district[state,2]+reffect.election[election,2]
  p1.old<-1/(1+exp(eta2.old)+exp(eta3.old))
  p2.old<-p1.old*exp(eta2.old)
  p3.old<-p1.old*exp(eta3.old)
 l.old<-log(p1.old)*(y==1)+log(p2.old)*(y==2)+log(p3.old)*(y==3)
 eta2.new<-X%*%b2+reffect.district.new[state,1]+reffect.election[election,1]
 eta3.new<-X%*%b3+reffect.district.new[state,2]+reffect.election[election,2]
  p1.new<-1/(1+exp(eta2.new)+exp(eta3.new))
  p2.new<-p1.new*exp(eta2.new)
  p3.new<-p1.new*exp(eta3.new)
 l.new<-log(p1.new)*(y==1)+log(p2.new)*(y==2)+log(p3.new)*(y==3)
 ratio<-(t(Matrix.District)%*%(l.new-l.old)) + 
	dmvnorm(reffect.district.new,rep(0,2),sigma.district, log=TRUE)-
	dmvnorm(reffect.district,rep(0,2),sigma.district, log=TRUE)
 u<-(log(runif(num.districts))<  ratio)*1
  reffect.district[!is.na(u) & u==1,]<-reffect.district.new[!is.na(u) & u==1,]
   reffect.district[,1]<-reffect.district[,1]-mean(reffect.district[,1])
   reffect.district[,2]<-reffect.district[,2]-mean(reffect.district[,2])

 ########################
 # Update sigma.district #
 #########################
sigma.district<-riwish(num.districts+3, .01*diag(2)+crossprod(reffect.district))
    

###################################
# Update Election reffects using MH #
 ###################################
reffect.election.new<-reffect.election+rmvt(num.elections,sigma=0.001*diag(2),3) 
  eta2.old<-X%*%b2+reffect.district[state,1]+reffect.election[election,1]
  eta3.old<-X%*%b3+reffect.district[state,2]+reffect.election[election,2]
  p1.old<-1/(1+exp(eta2.old)+exp(eta3.old))
  p2.old<-p1.old*exp(eta2.old)
  p3.old<-p1.old*exp(eta3.old)
 l.old<-log(p1.old)*(y==1)+log(p2.old)*(y==2)+log(p3.old)*(y==3)
 eta2.new<-X%*%b2+reffect.district[state,1]+reffect.election.new[election,1]
 eta3.new<-X%*%b3+reffect.district[state,2]+reffect.election.new[election,2]
  p1.new<-1/(1+exp(eta2.new)+exp(eta3.new))
  p2.new<-p1.new*exp(eta2.new)
  p3.new<-p1.new*exp(eta3.new)
 l.new<-log(p1.new)*(y==1)+log(p2.new)*(y==2)+log(p3.new)*(y==3)
 ratio<-(t(Matrix.Election)%*%(l.new-l.old)) + 
	dmvnorm(reffect.election.new,rep(0,2),sigma.election, log=TRUE)-
	dmvnorm(reffect.election,rep(0,2),sigma.election, log=TRUE)
 u<-(log(runif(num.elections))<  ratio)*1
  reffect.election[!is.na(u) & u==1,]<-reffect.election.new[!is.na(u) & u==1,]
   reffect.election[,1]<-reffect.election[,1]-mean(reffect.election[,1])
   reffect.election[,2]<-reffect.election[,2]-mean(reffect.election[,2])

 ##################
 # Update sigma.election #
 ##################
sigma.election<-riwish(num.elections+3, diag(2)+crossprod(reffect.election))

 
# Store Results
  if (i> burn & i%%thin==0) {
    j<-(i-burn)/thin
   Beta2[j,,ch]<-b2
   Beta3[j,,ch]<-b3
   Sigma.district[j,,ch]<-c(sigma.district)
   Sigma.election[j,,ch]<-c(sigma.election)
   Reffect.election[j,,ch]<-as.vector(reffect.election)
   Reffect.district[j,,ch]<-as.vector(reffect.district)
}

if (i%%500000==0)  print(c(paste("Chain=",ch), paste("Iteration=",i)))

# Backup saving (if something goes wrong)
if (i==nsim) {

Results<-list(Beta2=Beta2[,,ch], Beta3=Beta3[,,ch], Sigma.district=Sigma.district[,,ch], 
			Sigma.election=Sigma.election[,,ch], 
			Reffect.election=Reffect.election[,,ch],
			Reffect.district=Reffect.district[,,ch],
			vars=colnames(Beta2))
save(Results, file=paste("Individual-level model, chain=", ch))
}

}

} 
# End of MCMC algorithm


# Convergence checks
BETA.2<-mcmc.list(as.mcmc(Beta2[,,1]), as.mcmc(Beta2[,,2]), as.mcmc(Beta2[,,3]))
BETA.3<-mcmc.list(as.mcmc(Beta3[,,1]), as.mcmc(Beta3[,,2]), as.mcmc(Beta3[,,3]))
SIGMA.DISTRICT<-mcmc.list(as.mcmc(Sigma.district[,c(1,2,4),1]), as.mcmc(Sigma.district[,c(1,2,4),2]), as.mcmc(Sigma.district[,c(1,2,4),3]))
SIGMA.ELECTION<-mcmc.list(as.mcmc(Sigma.election[,c(1,2,4),1]), as.mcmc(Sigma.election[,c(1,2,4),2]), as.mcmc(Sigma.election[,c(1,2,4),3]))
gelman.diag(BETA.2)[[2]]<1.2
gelman.diag(BETA.3)[[2]]<1.2
gelman.diag(SIGMA.DISTRICT)[[2]]<1.2
gelman.diag(SIGMA.ELECTION)[[2]]<1.2

# Posterior summaries based on 1,000 convergent samples from each chain
Beta.2<-rbind(Beta2[9001:10000,,1], Beta2[9001:10000,,2], Beta2[9001:10000,,3])
Beta.3<-rbind(Beta3[9001:10000,,1], Beta3[9001:10000,,2], Beta3[9001:10000,,3])
Sigma.District<-rbind(Sigma.district[9001:10000,,1], Sigma.district[9001:10000,,2], Sigma.district[9001:10000,,3])
Sigma.Election<-rbind(Sigma.election[9001:10000,,1], Sigma.election[9001:10000,,2], Sigma.election[9001:10000,,3])
Reffect.Election<-rbind(Reffect.election[9001:10000,,1],Reffect.election[9001:10000,,2],Reffect.election[9001:10000,,3])
Reffect.District<-rbind(Reffect.district[9001:10000,,1],Reffect.district[9001:10000,,2],Reffect.district[9001:10000,,3])

Estimates.Table1<-list(Beta.2=Beta.2, Beta.3=Beta.3, Sigma.District=Sigma.District,
			Sigma.Election=Sigma.Election, Reffect.Election=Reffect.Election,
			Reffect.District=Reffect.District)
save(Estimates.Table1, file="Estimates.Table1.raw")


Coefficients.Absenteeism<-cbind(colMeans(Beta.3), HPDinterval(as.mcmc(Beta.3), p=0.9))
colnames(Coefficients.Absenteeism)<-c("Posterior mean", "90% CI - Low", "90% CI - High")

Coefficients.Invalid<-cbind(colMeans(Beta.2), HPDinterval(as.mcmc(Beta.2), p=0.9))
colnames(Coefficients.Invalid)<-c("Posterior mean", "90% CI - Low", "90% CI - High")

Table1.Col1<-rnd(Coefficients.Absenteeism[c(1,3,4,2,5,6,7:17),],2)
Table1.Col2<-rnd(Coefficients.Invalid[c(1,3,4,2,5,6,7:17),],2)

# Percent correctly predicted
mean.predicted.probs<-cbind(1/(1+
exp(X%*%apply(Beta.3,2,mean)+matrix(apply(Reffect.Election,2,mean),ncol=2)[election,2]+matrix(apply(Reffect.District,2,mean),ncol=2)[state,2])+
exp(X%*%apply(Beta.2,2,mean)+matrix(apply(Reffect.Election,2,mean),ncol=2)[election,1]+matrix(apply(Reffect.District,2,mean),ncol=2)[state,1])),
exp(X%*%apply(Beta.2,2,mean)+matrix(apply(Reffect.Election,2,mean),ncol=2)[election,1]+matrix(apply(Reffect.District,2,mean),ncol=2)[state,1])/(1+
exp(X%*%apply(Beta.3,2,mean)+matrix(apply(Reffect.Election,2,mean),ncol=2)[election,2]+matrix(apply(Reffect.District,2,mean),ncol=2)[state,2])+
exp(X%*%apply(Beta.2,2,mean)+matrix(apply(Reffect.Election,2,mean),ncol=2)[election,1]+matrix(apply(Reffect.District,2,mean),ncol=2)[state,1])) , 
exp(X%*%apply(Beta.3,2,mean)+matrix(apply(Reffect.Election,2,mean),ncol=2)[election,2]+matrix(apply(Reffect.District,2,mean),ncol=2)[state,2])/(1+
exp(X%*%apply(Beta.3,2,mean)+matrix(apply(Reffect.Election,2,mean),ncol=2)[election,2]+matrix(apply(Reffect.District,2,mean),ncol=2)[state,2])+
exp(X%*%apply(Beta.2,2,mean)+matrix(apply(Reffect.Election,2,mean),ncol=2)[election,1]+matrix(apply(Reffect.District,2,mean),ncol=2)[state,1])))
max.predicted.probs<-apply(mean.predicted.probs,1,max)
b<-sapply(1:nrow(X), function(i) which(mean.predicted.probs[i,]==max.predicted.probs[i]))
c<-sapply(1:nrow(X), function(i) which(Y[i,]==1))

# Chi-squared goodness of fit test 
# H0: Data is consistent with the fitted model
observed<-sapply(1:K, function(k) sum(Y[,k]))
expected<-sapply(1:K, function(k) sum(mean.predicted.probs[,k]))/nrow(X)

########################### ESTIMATES TABLE 1, COLS 1 & 2 ############### 
												#
#Column 1											#
Table1.Col1											#
												#
#Column 2											#
Table1.Col2											#
												#
# Percent correctly predicted								#
sum(b==c)/length(b)*100									#
												#
#Chi-squared test										#
chisq.test(observed, p=expected)$p.value  							#
												#
# Obs												#
nrow(Y)											#
#########################################################################

#Backup table in .csv
write.csv(Table1.Col1, file="Table1.Col1.csv")
write.csv(Table1.Col2, file="Table1.Col2.csv")

################################# FIGURE 3 #######################################################################################################################################################################################################
####### Expected change in the probabilities of non-voting associated with a  #####################################################################################################################################################################
####### change in the correlates of Cc, Cp, E, pB, and the efficiency of TREs #####################################################################################################################################################################

# Need to compute "marginal effects" first
dummies<-c(2,3,4,9,14,15,16)
categorical<-c(5,8,10)
continuous<-c(6,7,11,12,13,17)

Marginal.prob<-matrix(NA, nrow=(length(dummies)+length(categorical)+length(continuous))*2, ncol=5)
colnames(Marginal.prob)<-c("Dep.Variable", "Covariate", "Predictive Comparison",  "90% CI - Lower", "90% CI - Higher")
Marginal.prob[,1]<-c(rep("Absenteeism",(length(dummies)+length(categorical)+length(continuous))) ,rep("Invalid Voting",(length(dummies)+length(categorical)+length(continuous))))
Marginal.prob[,2]<-rep(colnames(X)[sort(c(dummies, categorical,continuous))],2)

for (j in dummies) {
X.0<-X
X.0[,j]<-0
X.1<-X
X.1[,j]<-1

Absenteeism<-sapply(1:dim(Estimates.Table1$Beta.2)[1], function(s) mean(
exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])/(1+
exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])) -
exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])/(1+
exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])))
)

Invalid<-sapply(1:dim(Estimates.Table1$Beta.3)[1], function(s) mean(
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])/(1+
exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])) -
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])/(1+
exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])))
)

Marginal.prob[grep(colnames(X)[j],Marginal.prob[,2]),3]<-c(mean(Absenteeism), mean(Invalid))
Marginal.prob[grep(colnames(X)[j],Marginal.prob[,2]),4:5]<-matrix(c(HPDinterval(as.mcmc(Absenteeism), p=0.9), HPDinterval(as.mcmc(Invalid), p=0.9)),ncol=2,byrow=TRUE)
}

for (j in categorical) {
X.0<-X
X.1<-X
X.1[,j]<-X.0[,j]+1

Absenteeism<-sapply(1:dim(Estimates.Table1$Beta.2)[1], function(s) mean(
exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])/(1+
exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])) -
exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])/(1+
exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])))
)

Invalid<-sapply(1:dim(Estimates.Table1$Beta.3)[1], function(s) mean(
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])/(1+
exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])) -
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])/(1+
exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])))
)

Marginal.prob[grep(colnames(X)[j],Marginal.prob[,2]),3]<-c(mean(Absenteeism), mean(Invalid))
Marginal.prob[grep(colnames(X)[j],Marginal.prob[,2]),4:5]<-matrix(c(HPDinterval(as.mcmc(Absenteeism), p=0.9), HPDinterval(as.mcmc(Invalid), p=0.9)),ncol=2,byrow=TRUE)
}

for (j in continuous) {
X.0<-X
X.1<-X
X.1[,j]<-X.0[,j]+sd(X[,j])

Absenteeism<-sapply(1:dim(Estimates.Table1$Beta.2)[1], function(s) mean(
exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])/(1+
exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])) -
exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])/(1+
exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])))
)

Invalid<-sapply(1:dim(Estimates.Table1$Beta.3)[1], function(s) mean(
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])/(1+
exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])) -
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])/(1+
exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])))
)

Marginal.prob[grep(colnames(X)[j],Marginal.prob[,2]),3]<-c(mean(Absenteeism), mean(Invalid))
Marginal.prob[grep(colnames(X)[j],Marginal.prob[,2]),4:5]<-matrix(c(HPDinterval(as.mcmc(Absenteeism), p=0.9), HPDinterval(as.mcmc(Invalid), p=0.9)),ncol=2,byrow=TRUE)
}

Marginal.prob<-data.frame(Marginal.prob)
for (k in 3:ncol(Marginal.prob)) {
Marginal.prob[,k]<-as.numeric(levels(Marginal.prob[,k]))[Marginal.prob[,k]]
}
names(Marginal.prob)<-c("Dep.Variable", "Covariate", "Mean","Lo.CI", "Hi.CI")
Marginal.prob[,3:5]<-100*Marginal.prob[,3:5]


Marginal.prob$Covariate<- factor(Marginal.prob$Covariate, 
	levels = rev(c("University", "Secondary", "Income",
		"Pol.Knowledge", "Urban", "Candidates",
		"Political Inefficacy", "No Political Representation",
		"Dissatisfaction w. Democracy", "Growth", "Inflation",
		"Clearance Rate", "Illiterates", "Young", "Seniors",
		"Competitiveness")),
	labels=rev(c("Education: University", "Education: Secondary",
		"Income", "Political Knowledge", 
		"Urban", "Candidates", "Political Inefficacy",
		"No Political Representation", "Dissatisfaction with Democracy",
		"Growth", "Inflation", "Clearance Rate",
		"Illiterates", "Young", "Seniors",
		"Competitiveness")))

Figure.3<-subset(Marginal.prob, Covariate!="Illiterates"  & Covariate!="Young" & 
			Covariate!="Seniors")
p.3<-ggplot(Figure.3,
 aes(y=Covariate, x=Mean)) + 
    geom_errorbarh(aes(xmin=Lo.CI, 
	xmax=Hi.CI, height=.2), colour="black") + facet_wrap(~Dep.Variable)+
    geom_point(size=3)+
theme_bw()+xlab("Marginal effect 
(in percentage points)")+ylab("")+
geom_hline(yintercept=0, linetype="longdash", colour="gray40")+
theme(strip.text.x = element_text(size = 12))+
geom_vline(xintercept=0, linetype="longdash", colour="gray40")

## Figure 3 - display #
	p.3

# Save figure 3
  jpeg("Figure3.jpg")
   p.3
   dev.off()
 
################################# FIGURE 4 #######################################################################################################################################################################################################
####### Expected change in the probabilities of non-voting associated with a  #####################################################################################################################################################################
####### change in the covariates capturing exemptions to compulsory voting #####################################################################################################################################################################
Figure.4<-subset(Marginal.prob, Covariate=="Illiterates"  | Covariate=="Young" | 
			Covariate=="Seniors")

p.4<-ggplot(Figure.4,
 aes(y=Covariate, x=Mean)) + 
    geom_errorbarh(aes(xmin=Lo.CI, 
	xmax=Hi.CI, height=.2), colour="black") + facet_wrap(~Dep.Variable)+
    geom_point(size=3)+
theme_bw()+xlab("Marginal effect 
(in percentage points)")+ylab("")+
geom_hline(yintercept=0, linetype="longdash", colour="gray40")+
theme(strip.text.x = element_text(size = 12))+
scale_x_continuous(breaks=c(-15,0, 15, 30, 45))+
geom_vline(xintercept=0, linetype="longdash", colour="gray40")


## Figure 4 - display #
	p.4

# Save figure 4
  jpeg("Figure4.jpg")
   p.4
   dev.off()

################################## TABLE 2 ######################################################################################################################################################################################################
################### MCMC (WinBUGS) estimates for the district-level model ###############################################################################################################################################################################
# Need to run WinBUGS model first

################ Columns 1 and 2 (baseline model) ######
df<-read.csv("district.1986.2014.csv")

# Create new vars
log_invalid<-log(df$Invalid/df$Valid)
log_absenteeism<-log(df$Absenteeism/df$Valid)
Y<-cbind(log_invalid, log_absenteeism)
n<-nrow(Y)
df$state<-as.numeric(factor(df$State), levels=unique(df$State))
df$election<-as.numeric(factor(df$Year), levels=unique(df$Year))
num.districts<-length(unique(df$state))
num.elections<-length(unique(df$election))

k<-2   # Number of outcome variables

# Priors 
R<-diag(2)
G0<-diag(2)
H0<-diag(2)
B<-diag(18)
b0=rep(0,18)
C<-diag(6)
c0=rep(0,6)


Compositional.Cols12.data<-list("num.districts"=num.districts,"num.elections"=num.elections,
	"n"=n, "B"=B, "C"=C, "b0"=b0, "c0"=c0,
 	"State"=df$state, "G0"=G0, "H0"=H0, "R"=R,
	 "Y"=Y, 
     "Urban"=df$Urban,
	"Income"=(df$Income/5),
	"Candidates"=(df$Candidates-mean(df$Candidates))/(2*sd(df$Candidates)),
	"Illiterate"=df$Illiterates, 
	 "Young"=df$Young, 
		"Old"=df$Seniors, 
	  "TRE"=df$TRE,
	 "Evote"=df$Electronic.Vote, 
	 "Competitiveness"=(df$Competitiveness-mean(df$Competitiveness))/(2*sd(df$Competitiveness)),
      "Election"=df$election,  
	"Growth"=(df$Growth-mean(df$Growth))/(2*sd(df$Growth)),
	"Inflation"=(log(df$Inflation)-mean(log(df$Inflation)))/(2*sd(log(df$Inflation)))
	  
)


set.seed(112215)  
Compositional.Cols12.inits<-function() {
list(T=diag(2), G=diag(2), H=diag(2), 
b=rnorm(18), c=rnorm(6),  
beta=matrix(rnorm(num.districts*2),num.districts,2),
alpha=matrix(rnorm(num.elections*2), num.elections,2),
k.e=round(runif(1,1,5)),
w.e=rgamma(n,1,1))} 

Compositional.Cols12.parameters<-c( "b", "c", "alpha", "beta", "w.e", "nu.e", 
	"Sigma", "Sigma.election", "Sigma.district")

Compositional.Cols12<-bugs(Compositional.Cols12.data,
	Compositional.Cols12.inits,Compositional.Cols12.parameters,
	"compositional.19862014.bug", 
n.iter=15000, n.burnin=5000,  
# Set the appropriate path for bugs.directory
bugs.directory="//isad.isadroot.ex.ac.uk/UOE/User/Desktop/WinBUGS14")
attach.bugs(Compositional.Cols12, overwrite=TRUE)

# Check convergence 
sum(Compositional.Cols12[[10]][,8]<1.2)/nrow(Compositional.Cols12[[10]])>0.975  #Neelon's convergence criterion

# Save estimates for tables in Online Appendix
save(Compositional.Cols12, file="Compositional.Cols12.raw")

#### Obtain posterior summaries for parameter estimates 
X<-cbind(1,	df$Urban,(df$Income/5),
		(df$Candidates-mean(df$Candidates))/(2*sd(df$Candidates)),
	df$Illiterates,df$Young,df$Seniors,df$TRE,
	df$Electronic.Vote,
	(df$Competitiveness-mean(df$Competitiveness))/(2*sd(df$Competitiveness)),
	(df$Growth-mean(df$Growth))/(2*sd(df$Growth)),	
	(log(df$Inflation)-mean(log(df$Inflation)))/(2*sd(log(df$Inflation))))

colnames(X)<-c("Intercept", "Urban", "Income", "Candidates","Illiterates", "Young",
			"Seniors", "Clearance Rate", "Electronic Vote",
			"Competitiveness",
			"Growth", "Inflation")

Estimates.Table2.Cols12<-matrix(NA, nrow=(ncol(X)*2), ncol=5)
colnames(Estimates.Table2.Cols12)<-c("Dep.Variable", "Covariate", "Posterior Mean",
	"90% CI, Low", "90% CI, High")
Estimates.Table2.Cols12[,1]<-c(rep("Absenteeism", ncol(X)),rep("Invalid Voting",ncol(X)))
Estimates.Table2.Cols12[,2]<-rep(colnames(X),2)

Estimates.Table2.Cols12[Estimates.Table2.Cols12[,1]=="Absenteeism",3:ncol(Estimates.Table2.Cols12)]<-rbind(c(mean(c[,4]), HPDinterval(as.mcmc(c[,4]),p=0.9)),
			cbind(apply(b[,10:18], 2,mean),HPDinterval(as.mcmc(b[,10:18]), p=0.9)),
			cbind(apply(c[,5:6], 2,mean),HPDinterval(as.mcmc(c[,5:6]), p=0.9)))

Estimates.Table2.Cols12[Estimates.Table2.Cols12[,1]=="Invalid Voting",3:ncol(Estimates.Table2.Cols12)]<-rbind(c(mean(c[,1]), HPDinterval(as.mcmc(c[,1]),p=0.9)),
		cbind(apply(b[,1:9], 2,mean),HPDinterval(as.mcmc(b[,1:9]), p=0.9)),
			cbind(apply(c[,2:3], 2,mean),HPDinterval(as.mcmc(c[,2:3]), p=0.9)))

Estimates.Table2.Cols12.df<-data.frame(Estimates.Table2.Cols12)
names(Estimates.Table2.Cols12.df)<-colnames(Estimates.Table2.Cols12)

for (k in 3:ncol(Estimates.Table2.Cols12.df)){
Estimates.Table2.Cols12.df[,k]<-rnd(as.numeric(levels(Estimates.Table2.Cols12.df[,k]))[Estimates.Table2.Cols12.df[,k]],2)
}

Table2.Col1<-subset(Estimates.Table2.Cols12.df, Dep.Variable=="Absenteeism")
Table2.Col2<-subset(Estimates.Table2.Cols12.df, Dep.Variable=="Invalid Voting")

########################### ESTIMATES TABLE 2, COLS 1 & 2 ##############
											     #	
#Column 1										     #
Table2.Col1[c(1,3,2,4,9,11,12,8,5,6,7,10),]				     #
											     #
#Column 2										     #
Table2.Col2[c(1,3,2,4,9,11,12,8,5,6,7,10),]				     #
											     #
# Deviance Information Criterion						     #
rnd(Compositional.Cols12[[10]][dim(Compositional.Cols12[[10]])][1],2)  #
											     #
#Observations									     #
n											     #
#######################################################################

#Backup table in .csv
write.csv(Table2.Col1[c(1,3,2,4,9,11,12,8,5,6,7,10),], file="Table2.Col1.csv")
write.csv(Table2.Col2[c(1,3,2,4,9,11,12,8,5,6,7,10),], file="Table2.Col2.csv")

################ Columns 3 and 4 
################ Interaction between Illiterates & Electronic Voting 

k<-2 # Number of outcome variables
# Priors 
R<-diag(2)
G0<-diag(2)
H0<-diag(2)
B<-diag(20)
b0=rep(0,20)
C<-diag(6)
c0=rep(0,6)


Compositional.Cols34.data<-list("num.districts"=num.districts,"num.elections"=num.elections,
	"n"=n, "B"=B, "C"=C, "b0"=b0, "c0"=c0,
 	"State"=df$state, "G0"=G0, "H0"=H0, "R"=R,
	 "Y"=Y, 
     "Urban"=df$Urban,
	"Income"=df$Income/5,
	"Candidates"=(df$Candidates-mean(df$Candidates))/(2*sd(df$Candidates)),
	"Illiterate"=df$Illiterates, 
	"Interaction.Illiterates"=df$Illiterates*df$Electronic.Vote,
	 "Young"=df$Young, 
		"Old"=df$Seniors, 
	  "TRE"=df$TRE,
	 "Evote"=df$Electronic.Vote,
	 "Competitiveness"=(df$Competitiveness-mean(df$Competitiveness))/(2*sd(df$Competitiveness)),
      "Election"=df$election,  
	"Growth"=(df$Growth-mean(df$Growth))/(2*sd(df$Growth)),
	"Inflation"=(log(df$Inflation)-mean(log(df$Inflation)))/(2*sd(log(df$Inflation)))
	  
)

set.seed(112215)  
Compositional.Cols34.inits<-function() {
list(T=diag(2), G=diag(2), H=diag(2), 
b=rnorm(20), c=rnorm(6),  
beta=matrix(rnorm(num.districts*2),num.districts,2),
alpha=matrix(rnorm(num.elections*2), num.elections,2),
k.e=round(runif(1,1,5)),
w.e=rgamma(n,1,1))} 


Compositional.Cols34.parameters<-c( "b", "c", "alpha", "beta", "w.e", "nu.e", 
	"Sigma", "Sigma.election", "Sigma.district")

Compositional.Cols34<-bugs(Compositional.Cols34.data,
	Compositional.Cols34.inits,Compositional.Cols34.parameters,
	"compositional.19862014.interaction.bug", 
n.iter=15000, n.burnin=5000,  
# Set the appropriate path for bugs.directory
bugs.directory="//isad.isadroot.ex.ac.uk/UOE/User/Desktop/WinBUGS14")
attach.bugs(Compositional.Cols34, overwrite=TRUE)

# Check convergence 
sum(Compositional.Cols34[[10]][,8]<1.2)/nrow(Compositional.Cols34[[10]])>0.975 #Neelon's criterion

# Save estimates for tables in Online Appendix
save(Compositional.Cols34, file="Compositional.Cols34.raw")

#### Obtain posterior summaries for parameter estimates 
X<-cbind(1,	df$Urban,df$Income/5,(df$Candidates-mean(df$Candidates))/(2*sd(df$Candidates)),
	df$Illiterates,df$Young,df$Seniors,df$TRE,
	df$Electronic.Vote,
	(df$Competitiveness-mean(df$Competitiveness))/(2*sd(df$Competitiveness)),
	df$Illiterates*df$Electronic.Vote,
	(df$Growth-mean(df$Growth))/(2*sd(df$Growth)),	
	(log(df$Inflation)-mean(log(df$Inflation)))/(2*sd(log(df$Inflation))))

colnames(X)<-c("Intercept", "Urban", "Income", "Candidates","Illiterates", "Young",
			"Seniors", "Clearance Rate", "Electronic Vote",
			"Competitiveness",
			"Interaction Illiterates-EVote",
			"Growth", "Inflation")

# Posterior Summaries
Estimates.Table2.Cols34<-matrix(NA, nrow=(ncol(X))*2, ncol=5)
colnames(Estimates.Table2.Cols34)<-c("Dep.Variable", "Covariate", "Posterior Mean",
	"90% CI, Low", "90% CI, High")
Estimates.Table2.Cols34[,1]<-c(rep("Absenteeism", (ncol(X))),rep("Invalid Voting",(ncol(X))))
Estimates.Table2.Cols34[,2]<-rep(colnames(X),2)

Estimates.Table2.Cols34[Estimates.Table2.Cols34[,1]=="Absenteeism",3:ncol(Estimates.Table2.Cols34)]<-rbind(c(mean(c[,4]),HPDinterval(as.mcmc(c[,4]), p=0.9)),
		cbind(apply(b[,c(10:18,20)], 2,mean),HPDinterval(as.mcmc(b[,c(10:18,20)]), p=0.9)),
			cbind(apply(c[,5:6], 2,mean),HPDinterval(as.mcmc(c[,5:6]), p=0.9)))

Estimates.Table2.Cols34[Estimates.Table2.Cols34[,1]=="Invalid Voting",3:ncol(Estimates.Table2.Cols34)]<-rbind(c(mean(c[,1]), HPDinterval(as.mcmc(c[,1]),p=0.9)),
			cbind(apply(b[,c(1:9,19)], 2,mean),HPDinterval(as.mcmc(b[,c(1:9,19)]), p=0.9)),
			cbind(apply(c[,2:3], 2,mean),HPDinterval(as.mcmc(c[,2:3]), p=0.9)))

Estimates.Table2.Cols34.df<-data.frame(Estimates.Table2.Cols34)
names(Estimates.Table2.Cols34.df)<-colnames(Estimates.Table2.Cols34)

for (k in 3:ncol(Estimates.Table2.Cols34.df)){
Estimates.Table2.Cols34.df[,k]<-rnd(as.numeric(levels(Estimates.Table2.Cols34.df[,k]))[Estimates.Table2.Cols34.df[,k]],2)
}

Table2.Col3<-subset(Estimates.Table2.Cols34.df, Dep.Variable=="Absenteeism")
Table2.Col4<-subset(Estimates.Table2.Cols34.df, Dep.Variable=="Invalid Voting")

########################### ESTIMATES TABLE 2, COLS 3 & 4 #################
												  #
# Column 3											  # 
Table2.Col3[c(1,3,2,4,9,12,13,8,5,11,6,7,10),-1]			        #
												  #
# Column 4											  # 											  #
Table2.Col4[c(1,3,2,4,9,12,13,8,5,11,6,7,10),-1]				  #
												  #
# Deviance Information Criterion							  #
rnd(Compositional.Cols34[[10]][dim(Compositional.Cols34[[10]])][1],2)     #
											        #
#Observations										  #
n												  #
												  #	
###########################################################################

#Backup table in .csv
write.csv(Table2.Col3[c(1,3,2,4,9,12,13,8,5,11,6,7,10),], file="Table2.Col3.csv")
write.csv(Table2.Col4[c(1,3,2,4,9,12,13,8,5,11,6,7,10),], file="Table2.Col4.csv")

################################# FIGURE 5 #######################################################################################################################################################################################################
####### Joint Impact of socio-demographic, electoral and political protest

# First need to compute the impact of each set of variables
# on the probabilities of Absenteeism & Invalid Votig
df.individual<-read.csv("individual.2002.2010.csv")
state<-as.numeric(factor(df.individual$State), levels=unique(df.individual$State))
election<-as.numeric(factor(df.individual$Year), levels=unique(df.individual$Year))
num.districts<-length(unique(state))
num.elections<-length(unique(election))
Y<-cbind(df.individual$Valid, df.individual$Invalid,df.individual$Absenteeism)
y<-c(sapply(1:nrow(Y), function(i) which(Y[i,]==1)))
n<-length(y)          # Number of obs
K<-3			    # Number of outcomes
X<-cbind(1,   df.individual$Income, df.individual$University,
		df.individual$Secondary,
		df.individual$Political.knowledge, df.individual$Urban,
		(df.individual$Candidates-mean(df.individual$Candidates))/(2*sd(df.individual$Candidates)), 
		df.individual$Political.inefficacy, 
		 df.individual$No.representation, 
		df.individual$Dissatisfaction.democracy, 
		(df.individual$Growth-mean(df.individual$Growth))/(2*sd(df.individual$Growth)),	
		(log(df.individual$Inflation)-mean(log(df.individual$Inflation)))/(2*sd(log(df.individual$Inflation))),
		df.individual$TRE, 
		df.individual$Illiterates, 
		df.individual$Young,
		df.individual$Seniors, 
		(df.individual$Competitiveness-mean(df.individual$Competitiveness))/(2*sd(df.individual$Competitiveness)),
		(df.individual$Age-mean(df.individual$Age))/(2*sd(df.individual$Age)),
		 df.individual$Female, df.individual$Married,
		df.individual$Religion, df.individual$Union, 
		df.individual$PT, df.individual$PSDB, df.individual$PMDB
		)
	  
colnames(X)<-c("Intercept", "Income", "University", "Secondary",
		"Pol.Knowledge",
			"Urban", "Candidates", "Political Inefficacy", "No Political Representation",
		"Dissatisfaction w. Democracy", "Growth", "Inflation",
		"Clearance Rate", "Illiterates", "Young", "Seniors", 
		"Competitiveness", 
	"Age", "Female", "Married", "Religion", "Union",
	 "PT", "PSDB", "PMDB")

###
X.0<-c(1,mean(X[,2]),1,0,max(X[,5]),max(X[,6]),min(X[,7]),
	mean(X[,8]),0,mean(X[,10]),mean(X[,11]),mean(X[,12]),
	mean(X[,13]),0,0,0,mean(X[,17]),mean(X[,18]),0,1,1,0,0,0,0)
X.1<-c(1,mean(X[,2]),0,0,min(X[,5]),min(X[,6]),max(X[,7]),
	mean(X[,8]),0,mean(X[,10]),mean(X[,11]),mean(X[,12]),
	mean(X[,13]),0,0,0,mean(X[,17]),mean(X[,18]),0,1,1,0,0,0,0)

Absenteeism.individual.demographicelectoral.1<-(sapply(1:dim(Estimates.Table1$Beta.2)[1], function(s) mean(
exp(X.1%*%Estimates.Table1$Beta.3[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[2]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[2])/(1+
exp(X.1%*%Estimates.Table1$Beta.3[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[2]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[2])+
exp(X.1%*%Estimates.Table1$Beta.2[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[1]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[1])))
)
)

Absenteeism.individual.demographicelectoral.0<-(sapply(1:dim(Estimates.Table1$Beta.2)[1], function(s) mean(
exp(X.0%*%Estimates.Table1$Beta.3[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[2]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[2])/(1+
exp(X.0%*%Estimates.Table1$Beta.3[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[2]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[2])+
exp(X.0%*%Estimates.Table1$Beta.2[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[1]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[1])))
)
)

Invalid.individual.demographicelectoral.1<-(sapply(1:dim(Estimates.Table1$Beta.2)[1], function(s) mean(
exp(X.1%*%Estimates.Table1$Beta.2[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[1]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[1])/(1+
exp(X.1%*%Estimates.Table1$Beta.3[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[2]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[2])+
exp(X.1%*%Estimates.Table1$Beta.2[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[1]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[1])))
)
)

Invalid.individual.demographicelectoral.0<-(sapply(1:dim(Estimates.Table1$Beta.2)[1], function(s) mean(
exp(X.0%*%Estimates.Table1$Beta.2[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[1]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[1])/(1+
exp(X.0%*%Estimates.Table1$Beta.3[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[2]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[2])+
exp(X.0%*%Estimates.Table1$Beta.2[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[1]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[1])))
)
)

X.1<-c(1,mean(X[,2]),0,0,min(X[,5]),min(X[,6]),max(X[,7]),
	mean(X[,8]),1,max(X[,10]),mean(X[,11]),mean(X[,12]),
	mean(X[,13]),0,0,0,mean(X[,17]),mean(X[,18]),0,1,1,0,0,0,0)

Absenteeism.individual.addpolitical.1<-(sapply(1:dim(Estimates.Table1$Beta.2)[1], function(s) mean(
exp(X.1%*%Estimates.Table1$Beta.3[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[2]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[2])/(1+
exp(X.1%*%Estimates.Table1$Beta.3[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[2]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[2])+
exp(X.1%*%Estimates.Table1$Beta.2[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[1]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[1])))
)
)

Absenteeism.individual.addpolitical.0<-(sapply(1:dim(Estimates.Table1$Beta.2)[1], function(s) mean(
exp(X.0%*%Estimates.Table1$Beta.3[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[2]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[2])/(1+
exp(X.0%*%Estimates.Table1$Beta.3[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[2]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[2])+
exp(X.0%*%Estimates.Table1$Beta.2[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[1]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[1])))
)
)

Invalid.individual.addpolitical.1<-(sapply(1:dim(Estimates.Table1$Beta.2)[1], function(s) mean(
exp(X.1%*%Estimates.Table1$Beta.2[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[1]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[1])/(1+
exp(X.1%*%Estimates.Table1$Beta.3[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[2]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[2])+
exp(X.1%*%Estimates.Table1$Beta.2[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[1]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[1])))
)
)

Invalid.individual.addpolitical.0<-(sapply(1:dim(Estimates.Table1$Beta.2)[1], function(s) mean(
exp(X.0%*%Estimates.Table1$Beta.2[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[1]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[1])/(1+
exp(X.0%*%Estimates.Table1$Beta.3[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[2]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[2])+
exp(X.0%*%Estimates.Table1$Beta.2[s,]+apply(matrix(Estimates.Table1$Reffect.Election[s,],ncol=2),2,mean)[1]+apply(matrix(Estimates.Table1$Reffect.District[s,],ncol=2),2,mean)[1])))
)
)

dd<-rbind(cbind("Absenteeism",c("Minimal-abstention", "Impact of socio-demographic 
	& electoral variables",
		"Impact of socio-demographic, 
        electoral & political protest variables"), 
c(mean(Absenteeism.individual.demographicelectoral.0),
mean(Absenteeism.individual.demographicelectoral.1),
mean(Absenteeism.individual.addpolitical.1))),
cbind("Invalid Voting",c("Minimal-abstention", "Impact of socio-demographic 
	& electoral variables",
		"Impact of socio-demographic, 
        electoral & political protest variables"),
c(mean(Invalid.individual.demographicelectoral.0),
mean(Invalid.individual.demographicelectoral.1),
mean(Invalid.individual.addpolitical.1)))
)

dd<-data.frame(dd)
names(dd)<-c("Dep.Variable", "Scenario", "Probability")
for (k in 3:3) {
dd[,k]<-as.numeric(levels(dd[,k]))[dd[,k]]
}
dd$Scenario<- factor(dd$Scenario, 
	levels = rev(c("Minimal-abstention",
			"Impact of socio-demographic 
	& electoral variables", 
	"Impact of socio-demographic, 
        electoral & political protest variables")))

p.5<-ggplot(dd, aes(x=Scenario, y=Probability))+
	geom_bar(stat="identity", fill="gray60", width=0.1)+
	facet_grid(.~Dep.Variable)+theme_bw()+
	ylab("Probabilities of non-voting")+xlab("")+
theme(strip.text.x = element_text(size = 12, face="bold"))+
theme(axis.title.x=element_text(vjust=-0.5))+
scale_y_continuous(breaks = c(0, 0.10, 0.20, 0.30, 0.40, 0.50))+
coord_flip()


### Figure 5 - display #
	p.5

# Save figure 5
  jpeg("Figure5.jpg")
   p.5
   dev.off()

################################ ONLINE APPENDIX ###############################################################################################################################################################################################################################
################################## TABLE A.1 #####################################################################################################################################################################################################################################
df<-read.csv("district.1986.2014.csv")
unique.state<-unique(df$State)
unique.year<-unique(df$Year)

#### Absenteeism by state
Absenteeism.state<-c()
for (j in 1:length(unique.state)) {
Absenteeism.state<-rbind(Absenteeism.state, 
c(as.character(unique.state[j]),
mean(df$Absenteeism[df$State==unique.state[j]]),
sd(df$Absenteeism[df$State==unique.state[j]]),
range(df$Absenteeism[df$State==unique.state[j]]))
)
}
Absenteeism.state<-data.frame(Absenteeism.state)
names(Absenteeism.state)<-c("State", "Mean", "Std.Dev", "Min", "Max")
for (k in 2:ncol(Absenteeism.state)) {
Absenteeism.state[,k]<-rnd(as.numeric(levels(Absenteeism.state[,k]))[Absenteeism.state[,k]],2)
}

#### Invalid voting by state
Invalid.state<-c()
for (j in 1:length(unique.state)) {
Invalid.state<-rbind(Invalid.state, 
c(as.character(unique.state[j]),
mean(df$Invalid[df$State==unique.state[j]]),
sd(df$Invalid[df$State==unique.state[j]]),
range(df$Invalid[df$State==unique.state[j]]))
)
}
Invalid.state<-data.frame(Invalid.state)
names(Invalid.state)<-c("State", "Mean", "Std.Dev", "Min", "Max")
for (k in 2:ncol(Invalid.state)) {
Invalid.state[,k]<-rnd(as.numeric(levels(Invalid.state[,k]))[Invalid.state[,k]],2)
}

#### Absenteeism by year
Absenteeism.year<-c()
for (j in 1:length(unique.year)) {
Absenteeism.year<-rbind(Absenteeism.year, 
c(as.character(unique.year[j]),
mean(df$Absenteeism[df$Year==unique.year[j]]),
sd(df$Absenteeism[df$Year==unique.year[j]]),
range(df$Absenteeism[df$Year==unique.year[j]]))
)
}
Absenteeism.year<-data.frame(Absenteeism.year)
names(Absenteeism.year)<-c("State", "Mean", "Std.Dev", "Min", "Max")
for (k in 2:ncol(Absenteeism.year)) {
Absenteeism.year[,k]<-rnd(as.numeric(levels(Absenteeism.year[,k]))[Absenteeism.year[,k]],2)
}

#### Invalid voting by year
Invalid.year<-c()
for (j in 1:length(unique.year)) {
Invalid.year<-rbind(Invalid.year, 
c(as.character(unique.year[j]),
mean(df$Invalid[df$Year==unique.year[j]]),
sd(df$Invalid[df$Year==unique.year[j]]),
range(df$Invalid[df$Year==unique.year[j]]))
)
}
Invalid.year<-data.frame(Invalid.year)
names(Invalid.year)<-c("State", "Mean", "Std.Dev", "Min", "Max")
for (k in 2:ncol(Invalid.year)) {
Invalid.year[,k]<-rnd(as.numeric(levels(Invalid.year[,k]))[Invalid.year[,k]],2)
}

#################################################
## UPPER-LEFT PANEL OF TABLE A.1:   		#
## Summary stats for Absenteeism by state       #
	Absenteeism.state					#
								#
## UPPER-RIGHT PANEL OF TABLE A.1:  		#
## Summary stats for Invalid voting by state    #
	Invalid.state					#
								#
## LOWER-LEFT PANEL OF TABLE A.1:   		#
## Summary stats for Absenteeism by year        #
	Absenteeism.year					#
								#
## LOWER-RIGHT PANEL OF TABLE A.1:   		#
## Summary stats for Absenteeism by year        #
	Invalid.year					#
#################################################

#Backup table in .csv
write.csv(Absenteeism.state, file="TableA1.AbsenteeismbyState.csv")
write.csv(Invalid.state, file="TableA1.InvalidbyState.csv")
write.csv(Absenteeism.year, file="TableA1.AbsenteeismbyYear.csv")
write.csv(Invalid.year, file="TableA1.InvalidbyYear.csv")

############################### FIGURE A.1 #####################################################################################################################################################################################################################################
df<-read.csv("district.1986.2014.csv")
df.2<-subset(df, select=c("State", "Absenteeism", "Invalid"))
df.2[,2]<-df.2[,2]*100
df.2[,3]<-df.2[,3]*100

# Multiple plot function
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
  library(grid)

  # Make a list from the ... arguments and plotlist
  plots <- c(list(...), plotlist)

  numPlots = length(plots)

  # If layout is NULL, then use 'cols' to determine layout
  if (is.null(layout)) {
    # Make the panel
    # ncol: Number of columns of plots
    # nrow: Number of rows needed, calculated from # of cols
    layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
                    ncol = cols, nrow = ceiling(numPlots/cols))
  }

 if (numPlots==1) {
    print(plots[[1]])

  } else {
    # Set up the page
    grid.newpage()
    pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))

    # Make each plot, in the correct location
    for (i in 1:numPlots) {
      # Get the i,j matrix positions of the regions that contain this subplot
      matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))

      print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
                                      layout.pos.col = matchidx$col))
    }
  }
}

df.Absenteeism<-c()
df.Invalid<-c()

for (i in 1:length(unique(df.2$State))) {
df.Absenteeism<-rbind(df.Absenteeism, c(as.character(unique(df.2$State)[i]),
quantile(df.2$Absenteeism[df.2$State==unique(df.2$State)[i]], 0.25),
quantile(df.2$Absenteeism[df.2$State==unique(df.2$State)[i]], 0.75),
median(df.2$Absenteeism[df.2$State==unique(df.2$State)[i]]),
quantile(df.2$Absenteeism[df.2$State==unique(df.2$State)[i]], 0.1),
quantile(df.2$Absenteeism[df.2$State==unique(df.2$State)[i]], 0.9)
)
)

df.Invalid<-rbind(df.Invalid, c(as.character(unique(df.2$State)[i]),
quantile(df.2$Invalid[df.2$State==unique(df.2$State)[i]], 0.25),
quantile(df.2$Invalid[df.2$State==unique(df.2$State)[i]], 0.75),
median(df.2$Invalid[df.2$State==unique(df.2$State)[i]]),
quantile(df.2$Invalid[df.2$State==unique(df.2$State)[i]], 0.1),
quantile(df.2$Invalid[df.2$State==unique(df.2$State)[i]], 0.9)
)
)

}

df.Absenteeism<-data.frame(df.Absenteeism)
names(df.Absenteeism)<-c("State", "lower", "upper", "middle", "ymin", "ymax")
for (k in 2:ncol(df.Absenteeism)) {
df.Absenteeism[,k]<-as.numeric(levels(df.Absenteeism[,k]))[df.Absenteeism[,k]]
}
df.Invalid<-data.frame(df.Invalid)
names(df.Invalid)<-c("State", "lower", "upper", "middle", "ymin", "ymax")
for (k in 2:ncol(df.Invalid)) {
df.Invalid[,k]<-as.numeric(levels(df.Invalid[,k]))[df.Invalid[,k]]
}

p1.A1 <- ggplot(df.Absenteeism, aes(x=State, lower=lower, upper=upper, middle=middle, ymin=ymin, ymax=ymax)) + 
geom_boxplot(stat="identity", outlier.shape=NA, fill="gray80")+theme_bw()+
xlab("State")+ggtitle("Absenteeism")+
theme(axis.title.y=element_text(vjust=1.25),
axis.title.x=element_text(vjust=-.35))+
ylab("As % of the electorate")+geom_hline(aes(yintercept=mean(df.2$Absenteeism)),linetype="dashed")

p2.A1 <- ggplot(df.Invalid, aes(x=State, lower=lower, upper=upper, middle=middle, ymin=ymin, ymax=ymax)) + 
geom_boxplot(stat="identity", outlier.shape=NA, fill="gray80")+theme_bw()+
xlab("State")+ggtitle("Invalid Voting")+
theme(axis.title.y=element_text(vjust=1.25),
axis.title.x=element_text(vjust=-.35))+
ylab("As % of the electorate")+geom_hline(aes(yintercept=mean(df.2$Invalid)),linetype="dashed")

# Figure A.1 - display #
multiplot(p1.A1, p2.A1, cols=1)

# Save figure A.1
  jpeg("FigureA1.jpg")
  multiplot(p1.A1, p2.A1, cols=1)
  dev.off()


################################# FIGURE A.2 #####################################################################################################################################################################################################################################
## Aggregate-level data
df<-read.csv("district.1986.2014.csv")

# Individual-level 
df.individual<-read.csv("individual.2002.2010.csv")

Comparison.2002<-cbind(sort(unique(as.character(df.individual$State))),
	as.numeric(by(df.individual$Absenteeism[df.individual$Year==2002],
		df.individual$State[df.individual$Year==2002],mean)),
as.numeric(by(df$Absenteeism[df$Year==2002],df$State[df$Year==2002],mean)),
as.numeric(by(df.individual$Invalid[df.individual$Year==2002],
		df.individual$State[df.individual$Year==2002],mean)),
as.numeric(by(df$Invalid[df$Year==2002],df$State[df$Year==2002],mean))
)

Comparison.2006<-cbind(sort(unique(as.character(df.individual$State))),
	as.numeric(by(df.individual$Absenteeism[df.individual$Year==2006],
		df.individual$State[df.individual$Year==2006],mean)),
as.numeric(by(df$Absenteeism[df$Year==2006],df$State[df$Year==2006],mean)),
as.numeric(by(df.individual$Invalid[df.individual$Year==2006],
		df.individual$State[df.individual$Year==2006],mean)),
as.numeric(by(df$Invalid[df$Year==2006],df$State[df$Year==2006],mean))
)

Comparison.2010<-cbind(sort(unique(as.character(df.individual$State))),
	as.numeric(by(df.individual$Absenteeism[df.individual$Year==2010],
		df.individual$State[df.individual$Year==2010],mean)),
as.numeric(by(df$Absenteeism[df$Year==2010],df$State[df$Year==2010],mean)),
as.numeric(by(df.individual$Invalid[df.individual$Year==2010],
		df.individual$State[df.individual$Year==2010],mean)),
as.numeric(by(df$Invalid[df$Year==2010],df$State[df$Year==2010],mean))
)

Comparison<-cbind(sort(unique(as.character(df.individual$State))),
	as.numeric(by(df.individual$Absenteeism,
		df.individual$State,mean)),
as.numeric(by(df$Absenteeism,df$State,mean)),
as.numeric(by(df.individual$Invalid,
		df.individual$State,mean)),
as.numeric(by(df$Invalid,df$State,mean))
)

# Figure A.2 - display #
par(mfrow=c(3,2))
par(mar=c(4.5,4.5,3,2))
plot(Comparison.2002[,2],Comparison.2002[,3], xlim=c(0,0.30),
	ylim=c(0,0.30), xlab=c("Individual-level data"),
	ylab=c("Official election returns"), main="Absenteeism, 2002")
abline(0,1)
plot(Comparison.2002[,4],Comparison.2002[,5], xlim=c(0,0.5),
	ylim=c(0,0.50),  xlab=c("Individual-level data"),
	ylab=c("Official election returns"), main="Invalid voting, 2002")
abline(0,1)
plot(Comparison.2006[,2],Comparison.2006[,3], xlim=c(0,0.40),
	ylim=c(0,0.40), xlab=c("Individual-level data"),
	ylab=c("Official election returns"), main="Absenteeism, 2006")
abline(0,1)
plot(Comparison.2006[,4],Comparison.2006[,5], xlim=c(0,1),
	ylim=c(0,1),  xlab=c("Individual-level data"),
	ylab=c("Official election returns"), main="Invalid voting, 2006")
abline(0,1)
plot(Comparison.2010[,2],Comparison.2010[,3], xlim=c(0,0.30),
	ylim=c(0,0.30), xlab=c("Individual-level data"),
	ylab=c("Official election returns"), main="Absenteeism, 2010")
abline(0,1)
plot(Comparison.2010[,4],Comparison.2010[,5], xlim=c(0,.6),
	ylim=c(0,.6),  xlab=c("Individual-level data"),
	ylab=c("Official election returns"), main="Invalid voting, 2010")
abline(0,1)

# Save figure A.2
  jpeg("FigureA2.jpg")
par(mfrow=c(3,2))
par(mar=c(4.5,4.5,3,2))
plot(Comparison.2002[,2],Comparison.2002[,3], xlim=c(0,0.30),
	ylim=c(0,0.30), xlab=c("Individual-level data"),
	ylab=c("Official election returns"), main="Absenteeism, 2002")
abline(0,1)
plot(Comparison.2002[,4],Comparison.2002[,5], xlim=c(0,0.5),
	ylim=c(0,0.50),  xlab=c("Individual-level data"),
	ylab=c("Official election returns"), main="Invalid voting, 2002")
abline(0,1)
plot(Comparison.2006[,2],Comparison.2006[,3], xlim=c(0,0.40),
	ylim=c(0,0.40), xlab=c("Individual-level data"),
	ylab=c("Official election returns"), main="Absenteeism, 2006")
abline(0,1)
plot(Comparison.2006[,4],Comparison.2006[,5], xlim=c(0,1),
	ylim=c(0,1),  xlab=c("Individual-level data"),
	ylab=c("Official election returns"), main="Invalid voting, 2006")
abline(0,1)
plot(Comparison.2010[,2],Comparison.2010[,3], xlim=c(0,0.30),
	ylim=c(0,0.30), xlab=c("Individual-level data"),
	ylab=c("Official election returns"), main="Absenteeism, 2010")
abline(0,1)
plot(Comparison.2010[,4],Comparison.2010[,5], xlim=c(0,.6),
	ylim=c(0,.6),  xlab=c("Individual-level data"),
	ylab=c("Official election returns"), main="Invalid voting, 2010")
abline(0,1)
  dev.off()

################################ TABLE A.2 #####################################################################################################################################################################################################################################
# Individual-level predictors
df.individual<-read.csv("individual.2002.2010.csv")
Individual.predictors<-rbind(c(rnd(mean(df.individual$Age),2), rnd(sd(df.individual$Age),2), rnd(range(df.individual$Age),2)),
c(rnd(mean(df.individual$Dissatisfaction.democracy),2), rnd(sd(df.individual$Dissatisfaction.democracy),2), rnd(range(df.individual$Dissatisfaction.democracy),2)),
c(rnd(mean(df.individual$University),2),rnd(sd(df.individual$University),2),rnd(range(df.individual$University),2)),
c(rnd(mean(df.individual$Secondary),2), rnd(sd(df.individual$Secondary),2), rnd(range(df.individual$Secondary),2)),
c(rnd(mean(df.individual$Female),2), rnd(sd(df.individual$Female),2), rnd(range(df.individual$Female),2)),
c(rnd(mean(df.individual$Illiterates),2),rnd(sd(df.individual$Illiterates),2),rnd(range(df.individual$Illiterates),2)),
c(rnd(mean(df.individual$Income),2),rnd(sd(df.individual$Income),2),rnd(range(df.individual$Income),2)),
c(rnd(mean(df.individual$Married),2),rnd(sd(df.individual$Married),2), rnd(range(df.individual$Married),2)),
c(rnd(mean(df.individual$No.representation),2),rnd(sd(df.individual$No.representation),2),rnd(range(df.individual$No.representation),2)),
c(rnd(mean(df.individual$PT),2),rnd(sd(df.individual$PT),2),rnd(range(df.individual$PT),2)),
c(rnd(mean(df.individual$PMDB),2), rnd(sd(df.individual$PMDB),2), rnd(range(df.individual$PMDB),2)),
c(rnd(mean(df.individual$PSDB),2), rnd(sd(df.individual$PSDB),2), rnd(range(df.individual$PSDB),2)),
c(rnd(mean(df.individual$Political.inefficacy),2),rnd(sd(df.individual$Political.inefficacy),2),rnd(range(df.individual$Political.inefficacy),2)),
c(rnd(mean(df.individual$Political.knowledge),2),rnd(sd(df.individual$Political.knowledge),2),rnd(range(df.individual$Political.knowledge),2)),
c(rnd(mean(df.individual$Religion),2), rnd(sd(df.individual$Religion),2), rnd(range(df.individual$Religion),2)),
c(rnd(mean(df.individual$Seniors),2),rnd(sd(df.individual$Seniors),2),rnd(range(df.individual$Seniors),2)),
c(rnd(mean(df.individual$Union),2), rnd(sd(df.individual$Union),2), rnd(range(df.individual$Union),2)),
c(rnd(mean(df.individual$Young),2), rnd(sd(df.individual$Young),2), rnd(range(df.individual$Young),2))
)
colnames(Individual.predictors)<-c("Mean", "Std. Dev", "Min", "Max")
rownames(Individual.predictors)<-c("(Log)Age", "Dissatisfaction with Democracy",
			"Education: University", "Education: Secondary", "Female",
			"Illiterates", "Income", "Married", "No Political Representation",
			"Partisanship:PT", "Partisanship:PMDB", "Partisanship:PSDB",
			 "Political Inefficacy", "Political Knowledge",
			"Religion", "Seniors", "Union", "Young")

## Aggregate-level predictors
df<-read.csv("district.1986.2014.csv")

# State-level predictors
State.predictors<-rbind(c(rnd(mean(df$Candidates),2),rnd(sd(df$Candidates),2),rnd(range(df$Candidates),2)),
c(rnd(mean(df$TRE),2),rnd(sd(df$TRE),2),rnd(range(df$TRE),2)),
c(rnd(mean(df$Competitiveness),2),rnd(sd(df$Competitiveness),2),rnd(range(df$Competitiveness),2)),
c(rnd(mean(df$Electronic.Vote),2),rnd(sd(df$Electronic.Vote),2),rnd(range(df$Electronic.Vote),2)),
c(rnd(mean(df$Illiterates),2), rnd(sd(df$Illiterates),2), rnd(range(df$Illiterates),2)),
c(rnd(mean(df$Income),2),rnd(sd(df$Income),2),rnd(range(df$Income),2)),
c(rnd(mean(df$Seniors),2),rnd(sd(df$Seniors),2), rnd(range(df$Seniors),2)),
c(rnd(mean(df$Urban),2), rnd(sd(df$Urban),2), rnd(range(df$Urban),2)),
c(rnd(mean(df$Young),2), rnd(sd(df$Young),2), rnd(range(df$Young),2)))
colnames(State.predictors)<-c("Mean", "Std. Dev", "Min", "Max")
rownames(State.predictors)<-c("Candidates", "Clearance Rate", "Competitiveness",
					"Electronic Vote", "Illiterates", "Income", "Seniors",
					"Urban", "Young")

# National-level predictors
National.predictors<-rbind(c(rnd(mean(by(df$Growth, df$Year, unique))*100,2),rnd(sd(by(df$Growth, df$Year, unique))*100,2),rnd(range(by(df$Growth, df$Year, unique))*100,2)),
c(rnd(mean(by(log(df$Inflation), df$Year, unique)),2),rnd(sd(by(log(df$Inflation), df$Year, unique)),2),rnd(range(by(log(df$Inflation), df$Year, unique)),2)))
colnames(National.predictors)<-c("Mean", "Std. Dev", "Min", "Max")
rownames(National.predictors)<-c("Growth(%)", "(Log)Inflation")


#####################################
## UPPER PANEL OF TABLE A.2:		#
  Individual.predictors			#
						#
## MIDDLE PANEL OF TABLE A.2:		#
  State.predictors			#
						#
## LOWER PANEL OF TABLE A.2:		#
  National.predictors			#
						#
#####################################

#Backup table in .csv
write.csv(Individual.predictors, file="TableA2.Individual.predictors.csv")
write.csv(State.predictors, file="TableA2.State.predictors.csv")
write.csv(National.predictors, file="TableA2.National.predictors.csv")

################################ FIGURE A.3 #####################################################################################################################################################################################################################################
dev.off()
load("Compositional.Cols12.raw")
attach.bugs(Compositional.Cols12, overwrite=TRUE)

# Figure A.3 - display #
split.screen(figs=c(2,1),erase=TRUE)
screen(1)
plot(density(nu.e), lwd=2.5, main="", xlab="", ylab="")
title(expression(paste("Posterior distribution of ", upsilon, " in the district-level model")), line=1)

screen(2)
hist(sapply(1:dim(w.e)[2], function(s) sum(w.e[,s]>=1)/length(w.e[,s])),
	main="", xlab="", ylab="")
title(expression(paste("Posterior probability Pr(", tau[j][t],">=1)")), line=2)


# Save figure A.3
 jpeg("FigureA3.jpg")
split.screen(figs=c(2,1),erase=TRUE)
screen(1)
plot(density(nu.e), lwd=2.5, main="", xlab="", ylab="")
title(expression(paste("Posterior distribution of ", upsilon, " in the district-level model")), line=1)

screen(2)
hist(sapply(1:dim(w.e)[2], function(s) sum(w.e[,s]>=1)/length(w.e[,s])),
	main="", xlab="", ylab="")
title(expression(paste("Posterior probability Pr(", tau[j][t],">=1)")), line=2)
  dev.off()

################################ FIGURE A.4 #####################################################################################################################################################################################################################################
dev.off()

# Need to compute residuals first
load("Compositional.Cols12.raw")
attach.bugs(Compositional.Cols12, overwrite=TRUE)
df<-read.csv("district.1986.2014.csv")

log_invalid<-log(df$Invalid/df$Valid)
log_absenteeism<-log(df$Absenteeism/df$Valid)
Y<-cbind(log_invalid, log_absenteeism)
n<-nrow(Y)
df$state<-as.numeric(factor(df$State), levels=unique(df$State))
df$election<-as.numeric(factor(df$Year), levels=unique(df$Year))
num.districts<-length(unique(df$state))
num.elections<-length(unique(df$election))

X<-cbind(1,	df$Urban,df$Income/5,
		(df$Candidates-mean(df$Candidates))/(2*sd(df$Candidates)),
	df$Illiterates,df$Young,df$Seniors,df$TRE,df$Electronic.Vote,
	(df$Competitiveness-mean(df$Competitiveness))/(2*sd(df$Competitiveness)),
	(df$Growth-mean(df$Growth))/(2*sd(df$Growth)),	
	(log(df$Inflation)-mean(log(df$Inflation)))/(2*sd(log(df$Inflation))))

colnames(X)<-c("Intercept", "Urban", "Income", "Candidates","Illiterates", "Young",
			"Seniors", "Clearance Rate", "Electronic Vote",
			"Competitiveness",
			"Growth", "Inflation")


# Univariate residuals (unweighted)
resid.I<-c()
resid.A<-c()
for (i in 1:nrow(df)) {
resid.I<-c(resid.I,
mean(sapply(1:dim(b)[1], function(s)
abs(Y[i,1]-X[i,]%*%c(c[s,1],as.vector(b[s,1:9]),as.vector(c[s,2:3])))/sqrt(Sigma[s,1,1])
))
)

resid.A<-c(resid.A,
mean(sapply(1:dim(b)[1], function(s)
abs(Y[i,2]-X[i,]%*%c(c[s,4],as.vector(b[s,10:18]),as.vector(c[s,5:6])))/sqrt(Sigma[s,2,2])
))
)
}

# Bivariate residuals (unweighted)
Resid<-c()
for (i in 1:nrow(df)) {
Resid<-c(Resid,
mean(sapply(1:dim(b)[1], function(s) cbind(Y[i,1]-X[i,]%*%c(c[s,1],as.vector(b[s,1:9]),as.vector(c[s,2:3])),
Y[i,2]-X[i,]%*%c(c[s,4],as.vector(b[s,10:18]),as.vector(c[s,5:6])))%*%
solve(Sigma[s,,])%*%
t(cbind(Y[i,1]-X[i,]%*%c(c[s,1],as.vector(b[s,1:9]),as.vector(c[s,2:3])),
Y[i,2]-X[i,]%*%c(c[s,4],as.vector(b[s,10:18]),as.vector(c[s,5:6]))))))
)
}

# Figure A.4 - display #
split.screen(figs=c(2,1),erase=TRUE)

split.screen(c(1, 2), screen = 1)
screen(3)
 plot(resid.A, col="gray", ylab="Unweighted residuals", 
xlab="Observations")
title(expression(paste(abs(epsilon[j][t]^A)," / ", sigma[epsilon]^A)),line=1) 
abline(h=3, lwd=2)

screen(4)
 plot(resid.I, col="gray", ylab="Unweighted residuals", 
xlab="Observations")
abline(h=3, lwd=2)
title(expression(paste(abs(epsilon[j][t]^I)," / ", sigma[epsilon]^I)),line=1) 

screen(2)
plot(Resid, col="gray",ylab="Unweighted residuals", 
xlab="Observations")
abline(h=qchisq((1-pnorm(-3)*2), df=2), lwd=2)
title(expression(paste(epsilon[j][t]^T, Sigma[epsilon]^-1,epsilon[j][t])),line=1 )


dev.off()

# Save figure A.4
 jpeg("FigureA4.jpg")
split.screen(figs=c(2,1),erase=TRUE)
split.screen(c(1, 2), screen = 1)
screen(3)
 plot(resid.A, col="gray", ylab="Unweighted residuals", 
xlab="Observations")
title(expression(paste(abs(epsilon[j][t]^A)," / ", sigma[epsilon]^A)),line=1) 
abline(h=3, lwd=2)

screen(4)
 plot(resid.I, col="gray", ylab="Unweighted residuals", 
xlab="Observations")
abline(h=3, lwd=2)
title(expression(paste(abs(epsilon[j][t]^I)," / ", sigma[epsilon]^I)),line=1) 

screen(2)
plot(Resid, col="gray",ylab="Unweighted residuals", 
xlab="Observations")
abline(h=qchisq((1-pnorm(-3)*2), df=2), lwd=2)
title(expression(paste(epsilon[j][t]^T, Sigma[epsilon]^-1,epsilon[j][t])),line=1 )
dev.off()

################################ FIGURE A.5 #####################################################################################################################################################################################################################################

# For Figure A.5, we need to estimate pooled OLS and HLM models
df<-read.csv("district.1986.2014.csv")

# OLS estimation
n<-length(df$Invalid)
B<-diag(18)
C<-diag(6)

PooledOLS.data<-list("n"=n, "B"=B, "C"=C,
 	 "Invalid"=df$Invalid,
	 "Absenteeism"=df$Absenteeism ,
     "Urban"=df$Urban,
	"Income"=df$Income/5,
	"Candidates"=(df$Candidates-mean(df$Candidates))/(2*sd(df$Candidates)),
	"Illiterate"=df$Illiterates, 
	 "Young"=df$Young, 
		"Old"=df$Seniors, 
	  "TRE"=df$TRE,
	 "Evote"=df$Electronic.Vote, 
	 "Competitiveness"=(df$Competitiveness-mean(df$Competitiveness))/(2*sd(df$Competitiveness)),
	"Growth"=(df$Growth-mean(df$Growth))/(2*sd(df$Growth)),
	"Inflation"=(log(df$Inflation)-mean(log(df$Inflation)))/(2*sd(log(df$Inflation)))
)

set.seed(112215)  
PooledOLS.inits<-function() {
list(b=rnorm(18), c=rnorm(6),
tau.invalid=runif(1,0,1),
tau.absenteeism=runif(1,0,1))} 

PooledOLS.parameters<-c( "b", "c",  "sigma.invalid", "sigma.absenteeism")

PooledOLS<-bugs(PooledOLS.data,
	PooledOLS.inits,PooledOLS.parameters,
	"PooledOLS.bug", 
n.iter=15000, n.burnin=5000,  
# Set the appropriate path for bugs.directory
bugs.directory="//isad.isadroot.ex.ac.uk/UOE/User/Desktop/WinBUGS14")
attach.bugs(PooledOLS, overwrite=TRUE)

# Check convergence
sum(PooledOLS[[10]][,8]<1.2)/nrow(PooledOLS[[10]])>0.975 #Neelon's criterion

# Save estimates for use in Figure A.5 and Tables A.3 & A.4
save(PooledOLS, file="PooledOLS.raw")

X<-cbind(1,	df$Urban,df$Income/5,(df$Candidates-mean(df$Candidates))/(2*sd(df$Candidates)),
	df$Illiterates,df$Young,df$Seniors,df$TRE,df$Electronic.Vote,
	(df$Competitiveness-mean(df$Competitiveness))/(2*sd(df$Competitiveness)),
	(df$Growth-mean(df$Growth))/(2*sd(df$Growth)),	
	(log(df$Inflation)-mean(log(df$Inflation)))/(2*sd(log(df$Inflation))))

colnames(X)<-c("Intercept", "Urban", "Income", "Candidates","Illiterates", "Young",
			"Seniors", "Clearance Rate", "Electronic Vote",
			"Competitiveness",
			"Growth", "Inflation")

# POSTERIOR SIMULATIONS
P.simulated.OLS<-array(NA,c(n,2,PooledOLS[[6]]))
for (i in 1:n) {
P.simulated.OLS[i,1,]<-sapply(1:PooledOLS[[6]], function(s) rnorm(1,
c(X[i,]%*%c(c[s,1],as.vector(b[s,1:9]),as.vector(c[s,2:3]))),
sigma.invalid[s]))

P.simulated.OLS[i,2,]<-sapply(1:PooledOLS[[6]], function(s) rnorm(1,
X[i,]%*%c(c[s,4],as.vector(b[s,10:18]),as.vector(c[s,5:6])),
sigma.absenteeism[s]))
}

SimulatedOLS<-data.frame(as.vector(P.simulated.OLS[,1,]), as.vector(P.simulated.OLS[,2,]))
names(SimulatedOLS)<-c("Invalid Voting", "Absenteeism")


# HLM estimation
df$state<-as.numeric(factor(df$State), levels=unique(df$State))
df$election<-as.numeric(factor(df$Year), levels=unique(df$Year))
num.districts<-length(unique(df$state))
num.elections<-length(unique(df$election))

B<-diag(18)
C<-diag(6)

Multilevelregression.data<-list("n"=n, "B"=B, "C"=C,
"State"=df$state,   "Election"=df$election,
 	 "Invalid"=df$Invalid,
"num.districts"=num.districts,"num.elections"=num.elections,
	 "Absenteeism"=df$Absenteeism ,
     "Urban"=df$Urban,
	"Income"=df$Income/5,
	"Candidates"=(df$Candidates-mean(df$Candidates))/(2*sd(df$Candidates)),
	"Illiterate"=df$Illiterates, 
	 "Young"=df$Young, 
		"Old"=df$Seniors, 
	  "TRE"=df$TRE,
	 "Evote"=df$Electronic.Vote, 
	 "Competitiveness"=(df$Competitiveness-mean(df$Competitiveness))/(2*sd(df$Competitiveness)),
	"Growth"=(df$Growth-mean(df$Growth))/(2*sd(df$Growth)),
	"Inflation"=(log(df$Inflation)-mean(log(df$Inflation)))/(2*sd(log(df$Inflation)))
)

set.seed(112215)  
Multilevelregression.inits<-function() {
list(b=rnorm(18), c=rnorm(6),
tau.invalid=runif(1,0,1),
tau.absenteeism=runif(1,0,1),
tau.district=c(1,1),
beta=matrix(rnorm(num.districts*2),num.districts,2),
alpha=matrix(rnorm(num.elections*2), num.elections,2),
tau.election=c(1,1))} 

Multilevelregression.parameters<-c( "b", "c",  "sigma.invalid", 
"sigma.absenteeism", "sigma.election", "sigma.district", "beta", "alpha")

Multilevelregression<-bugs(Multilevelregression.data,
	Multilevelregression.inits,Multilevelregression.parameters,
	"Multilevelregression.bug", 
n.iter=35000, n.burnin=5000,  
# Set the appropriate path for bugs.directory
bugs.directory="//isad.isadroot.ex.ac.uk/UOE/User/Desktop/WinBUGS14")
attach.bugs(Multilevelregression, overwrite=TRUE)

# Check convergence
sum(Multilevelregression[[10]][,8]<1.2)/nrow(Multilevelregression[[10]])>0.975 #Neelon's criterion

# Save estimates for use in Figure A.5 and Tables A.3 & A.4
save(Multilevelregression, file="Multilevelregression.raw")

### POSTERIOR SUMULATIONS
X<-cbind(1,	df$Urban,df$Income/5,(df$Candidates-mean(df$Candidates))/(2*sd(df$Candidates)),
	df$Illiterates,df$Young,df$Seniors,df$TRE,df$Electronic.Vote,
	(df$Competitiveness-mean(df$Competitiveness))/(2*sd(df$Competitiveness)),
	(df$Growth-mean(df$Growth))/(2*sd(df$Growth)),	
	(log(df$Inflation)-mean(log(df$Inflation)))/(2*sd(log(df$Inflation))))

colnames(X)<-c("Intercept", "Urban", "Income", "Candidates","Illiterates", "Young",
			"Seniors", "Clearance Rate", "Electronic Vote",
			"Competitiveness",
			"Growth", "Inflation")

P.simulated.HLM<-array(NA,c(n,2,Multilevelregression[[6]]))
elections.simul<-array(NA,c(num.elections, 2,Multilevelregression[[6]]))
districts.simul<-array(NA,c(num.districts, 2,Multilevelregression[[6]]))

for (i in 1:num.elections) {
elections.simul[i,1,]<-sapply(1:Multilevelregression[[6]], function(s) rnorm(1,c(0,0),sigma.election[s,1]))
elections.simul[i,2,]<-sapply(1:Multilevelregression[[6]], function(s) rnorm(1,c(0,0),sigma.election[s,2]))
}

for (i in 1:num.districts) {
districts.simul[i,1,]<-sapply(1:Multilevelregression[[6]], function(s) rnorm(1,c(0,0),sigma.district[s,1]))
districts.simul[i,2,]<-sapply(1:Multilevelregression[[6]], function(s) rnorm(1,c(0,0),sigma.district[s,2]))
}

for (i in 1:n) {
P.simulated.HLM[i,1,]<-sapply(1:Multilevelregression[[6]], function(s) rnorm(1,
c(X[i,]%*%c(c[s,1],as.vector(b[s,1:9]),as.vector(c[s,2:3])))+
elections.simul[df$election[i],1,s]+districts.simul[df$state[i],1,s],
sigma.invalid[s]))

P.simulated.HLM[i,2,]<-sapply(1:Multilevelregression[[6]], function(s) rnorm(1,
X[i,]%*%c(c[s,4],as.vector(b[s,10:18]),as.vector(c[s,5:6]))+
 elections.simul[df$election[i],2,s]+districts.simul[df$state[i],2,s],
sigma.absenteeism[s]))
}

SimulatedHLM<-data.frame(as.vector(P.simulated.HLM[,1,]), as.vector(P.simulated.HLM[,2,]))
names(SimulatedHLM)<-c("Invalid Voting", "Absenteeism")

# Figure A.5 - display #
split.screen(figs=c(2,1),erase=TRUE)
screen(1)
 plot(SimulatedOLS[,1], SimulatedOLS[,2],
	col="gray", ylab="Absenteeism (%)",  
xlab="Invalid Voting (%)")
title("Predictions from pooled OLS models",line=1) 
abline(h=0, lwd=2)
abline(v=0, lwd=2)

screen(2)
 plot(SimulatedHLM[,1], SimulatedHLM[,2],
	col="gray", ylab="Absenteeism (%)",  
xlab="Invalid Voting (%)")
title("Predictions from hierarchical linear models",line=1) 
abline(h=0, lwd=2)
abline(v=0, lwd=2)


dev.off()

# Save figure A.5
 jpeg("FigureA5.jpg")
split.screen(figs=c(2,1),erase=TRUE)
screen(1)
 plot(SimulatedOLS[,1], SimulatedOLS[,2],
	col="gray", ylab="Absenteeism (%)",  
xlab="Invalid Voting (%)")
title("Predictions from pooled OLS models",line=1) 
abline(h=0, lwd=2)
abline(v=0, lwd=2)

screen(2)
 plot(SimulatedHLM[,1], SimulatedHLM[,2],
	col="gray", ylab="Absenteeism (%)",  
xlab="Invalid Voting (%)")
title("Predictions from hierarchical linear models",line=1) 
abline(h=0, lwd=2)
abline(v=0, lwd=2)
dev.off()

################################ TABLE A.3 #####################################################################################################################################################################################################################################

df<-read.csv("district.1986.2014.csv")
X<-cbind(1,	df$Urban,df$Income/5,(df$Candidates-mean(df$Candidates))/(2*sd(df$Candidates)),
	df$Illiterates,df$Young,df$Seniors,df$TRE,df$Electronic.Vote,
	(df$Competitiveness-mean(df$Competitiveness))/(2*sd(df$Competitiveness)),
	(df$Growth-mean(df$Growth))/(2*sd(df$Growth)),	
	(log(df$Inflation)-mean(log(df$Inflation)))/(2*sd(log(df$Inflation))))

colnames(X)<-c("Intercept", "Urban", "Income", "Candidates","Illiterates", "Young",
			"Seniors", "Clearance Rate", "Electronic Vote",
			"Competitiveness",
			"Growth", "Inflation")

# Posterior Summaries OLS
load("PooledOLS.raw")
attach.bugs(PooledOLS, overwrite=TRUE)
Estimates.OLS<-matrix(NA, nrow=(ncol(X)*2), ncol=5)
colnames(Estimates.OLS)<-c("Dep.Variable", "Covariate", "Posterior Mean",
	"90% CI, Low", "90% CI")
Estimates.OLS[,1]<-c(rep("Absenteeism", ncol(X)),rep("Invalid Voting",ncol(X)))
Estimates.OLS[,2]<-rep(colnames(X),2)

Estimates.OLS[Estimates.OLS[,1]=="Absenteeism",3:ncol(Estimates.OLS)]<-rbind(c(mean(c[,4]), HPDinterval(as.mcmc(c[,4]),p=0.9)),
			cbind(apply(b[,10:18], 2,mean),HPDinterval(as.mcmc(b[,10:18]), p=0.9)),
			cbind(apply(c[,5:6], 2,mean),HPDinterval(as.mcmc(c[,5:6]), p=0.9)))

Estimates.OLS[Estimates.OLS[,1]=="Invalid Voting",3:ncol(Estimates.OLS)]<-rbind(c(mean(c[,1]), HPDinterval(as.mcmc(c[,1]),p=0.9)),
		cbind(apply(b[,1:9], 2,mean),HPDinterval(as.mcmc(b[,1:9]), p=0.9)),
			cbind(apply(c[,2:3], 2,mean),HPDinterval(as.mcmc(c[,2:3]), p=0.9)))

Estimates.OLS.df<-data.frame(Estimates.OLS)
names(Estimates.OLS.df)<-colnames(Estimates.OLS)
for (k in 3:ncol(Estimates.OLS.df)){
Estimates.OLS.df[,k]<-rnd(as.numeric(levels(Estimates.OLS.df[,k]))[Estimates.OLS.df[,k]],2)
}

# Posterior summaries HLM
load("Multilevelregression.raw")
attach.bugs(Multilevelregression, overwrite=TRUE)

Estimates.HLM<-matrix(NA, nrow=(ncol(X)*2), ncol=5)
colnames(Estimates.HLM)<-c("Dep.Variable", "Covariate", "Posterior Mean",
	"90% CI, Low", "90% CI")
Estimates.HLM[,1]<-c(rep("Absenteeism", ncol(X)),rep("Invalid Voting",ncol(X)))
Estimates.HLM[,2]<-rep(colnames(X),2)

Estimates.HLM[Estimates.HLM[,1]=="Absenteeism",3:ncol(Estimates.OLS)]<-rbind(c(mean(c[,4]), HPDinterval(as.mcmc(c[,4]),p=0.9)),
			cbind(apply(b[,10:18], 2,mean),HPDinterval(as.mcmc(b[,10:18]), p=0.9)),
			cbind(apply(c[,5:6], 2,mean),HPDinterval(as.mcmc(c[,5:6]), p=0.9)))

Estimates.HLM[Estimates.HLM[,1]=="Invalid Voting",3:ncol(Estimates.OLS)]<-rbind(c(mean(c[,1]), HPDinterval(as.mcmc(c[,1]),p=0.9)),
		cbind(apply(b[,1:9], 2,mean),HPDinterval(as.mcmc(b[,1:9]), p=0.9)),
			cbind(apply(c[,2:3], 2,mean),HPDinterval(as.mcmc(c[,2:3]), p=0.9)))

Estimates.HLM.df<-data.frame(Estimates.HLM)
names(Estimates.HLM.df)<-colnames(Estimates.HLM)

for (k in 3:ncol(Estimates.HLM.df)){
Estimates.HLM.df[,k]<-rnd(as.numeric(levels(Estimates.HLM.df[,k]))[Estimates.HLM.df[,k]],2)
}

Table.A3.Col1<-subset(Estimates.OLS.df,  Dep.Variable=="Absenteeism")
Table.A3.Col2<-subset(Estimates.OLS.df,  Dep.Variable=="Invalid Voting")
Table.A3.Col3<-subset(Estimates.HLM.df,  Dep.Variable=="Absenteeism")
Table.A3.Col4<-subset(Estimates.HLM.df,  Dep.Variable=="Invalid Voting")

######################## TABLE A.3 ####################
									#
# Column 1								#
Table.A3.Col1[c(1,3,2,4,9,8,5,6,7,11,12,10),] 		#
									#
# Column 2								#
Table.A3.Col2[c(1,3,2,4,9,8,5,6,7,11,12,10),]		#
									#
# Column 3								#
Table.A3.Col3[c(1,3,2,4,9,8,5,6,7,11,12,10),]		#
									#
# Column 4								#
Table.A3.Col4[c(1,3,2,4,9,8,5,6,7,11,12,10),]		#
									#
#######################################################

#Backup table in .csv
write.csv(Table.A3.Col1[c(1,3,2,4,9,8,5,6,7,11,12,10),], file="TableA3.Col1.csv")
write.csv(Table.A3.Col2[c(1,3,2,4,9,8,5,6,7,11,12,10),], file="TableA3.Col2.csv")
write.csv(Table.A3.Col3[c(1,3,2,4,9,8,5,6,7,11,12,10),], file="TableA3.Col3.csv")
write.csv(Table.A3.Col4[c(1,3,2,4,9,8,5,6,7,11,12,10),], file="TableA3.Col4.csv")


################################ TABLE A.4 #####################################################################################################################################################################################################################################
df<-read.csv("district.1986.2014.csv")

log_invalid<-log(df$Invalid/df$Valid)
log_absenteeism<-log(df$Absenteeism/df$Valid)
Y<-cbind(log_invalid, log_absenteeism)
n<-nrow(Y)
df$state<-as.numeric(factor(df$State), levels=unique(df$State))
df$election<-as.numeric(factor(df$Year), levels=unique(df$Year))
num.districts<-length(unique(df$state))
num.elections<-length(unique(df$election))

X<-cbind(1,	df$Urban,df$Income/5,(df$Candidates-mean(df$Candidates))/(2*sd(df$Candidates)),
	df$Illiterates,df$Young,df$Seniors,df$TRE,df$Electronic.Vote,
	(df$Competitiveness-mean(df$Competitiveness))/(2*sd(df$Competitiveness)),
	(df$Growth-mean(df$Growth))/(2*sd(df$Growth)),	
	(log(df$Inflation)-mean(log(df$Inflation)))/(2*sd(log(df$Inflation))))

colnames(X)<-c("Intercept", "Urban", "Income", "Candidates","Illiterates", "Young",
			"Seniors", "Clearance Rate", "Electronic Vote",
			"Competitiveness",
			"Growth", "Inflation")
n<-nrow(X)

# Need to calculate d1 and d2 first
## Compositional model
load("Compositional.Cols12.raw")
attach.bugs(Compositional.Cols12, overwrite=TRUE)

Y.simulated<-array(NA,c(n,2,Compositional.Cols12[[6]]))
Y.expected<-array(NA,c(n,2,Compositional.Cols12[[6]]))
P.simulated<-array(NA,c(n,2,Compositional.Cols12[[6]]))
P.expected<-array(NA,c(n,2,Compositional.Cols12[[6]]))
elections.simul<-array(NA,c(num.elections, 2,Compositional.Cols12[[6]]))
districts.simul<-array(NA,c(num.districts, 2,Compositional.Cols12[[6]]))

for (i in 1:num.elections) {
elections.simul[i,,]<-sapply(1:Compositional.Cols12[[6]], function(s) mvrnorm(1,c(0,0),Sigma.election[s,,]))
}
for (i in 1:num.districts) {
districts.simul[i,,]<-sapply(1:Compositional.Cols12[[6]], function(s) mvrnorm(1,c(0,0),Sigma.district[s,,]))
}
for (i in 1:n) {
Y.simulated[i,,]<-sapply(1:Compositional.Cols12[[6]], function(s) mvrnorm(1,
c(X[i,]%*%c(c[s,1],as.vector(b[s,1:9]),as.vector(c[s,2:3]))+
 elections.simul[df$election[i],1,s]+districts.simul[df$state[i],1,s],
X[i,]%*%c(c[s,4],as.vector(b[s,10:18]),as.vector(c[s,5:6]))+
 elections.simul[df$election[i],2,s]+districts.simul[df$state[i],2,s]),
(1/w.e[s,i])*Sigma[s,,]))

Y.expected[i,,]<-sapply(1:Compositional.Cols12[[6]], function(s) c(X[i,]%*%c(c[s,1],as.vector(b[s,1:9]),as.vector(c[s,2:3])),
X[i,]%*%c(c[s,4],as.vector(b[s,10:18]),as.vector(c[s,5:6]))))

P.simulated[i,,]<-cbind(exp(Y.simulated[i,1,])/(1+exp(Y.simulated[i,1,])+exp(Y.simulated[i,2,])),
exp(Y.simulated[i,2,])/(1+exp(Y.simulated[i,1,])+exp(Y.simulated[i,2,])))

P.expected[i,,]<-cbind(exp(Y.expected[i,1,])/(1+exp(Y.expected[i,1,])+exp(Y.expected[i,2,])),
exp(Y.expected[i,2,])/(1+exp(Y.expected[i,1,])+exp(Y.expected[i,2,])))
}

d1.compositional<-sapply(1:Compositional.Cols12[[6]],function(s) (P.expected[,1,s]-df$Invalid)^2+(P.expected[,2,s]-df$Absenteeism)^2)
d2.compositional<-sapply(1:Compositional.Cols12[[6]], function(s)(df$Invalid-P.simulated[,1,s])^2+(df$Absenteeism-P.simulated[,2,s])^2)

## HLM model
load("Multilevelregression.raw")
attach.bugs(Multilevelregression, overwrite=TRUE)

P.simulated<-array(NA,c(n,2,Multilevelregression[[6]]))
P.expected<-array(NA,c(n,2,Multilevelregression[[6]]))
elections.simul<-array(NA,c(num.elections, 2,Multilevelregression[[6]]))
districts.simul<-array(NA,c(num.districts, 2,Multilevelregression[[6]]))

for (i in 1:num.elections) {
elections.simul[i,1,]<-sapply(1:Multilevelregression[[6]], function(s) rnorm(1,c(0,0),sigma.election[s,1]))
elections.simul[i,2,]<-sapply(1:Multilevelregression[[6]], function(s) rnorm(1,c(0,0),sigma.election[s,2]))
}
for (i in 1:num.districts) {
districts.simul[i,1,]<-sapply(1:Multilevelregression[[6]], function(s) rnorm(1,c(0,0),sigma.district[s,1]))
districts.simul[i,2,]<-sapply(1:Multilevelregression[[6]], function(s) rnorm(1,c(0,0),sigma.district[s,2]))
}
for (i in 1:n) {
P.simulated[i,1,]<-sapply(1:Multilevelregression[[6]], function(s) rnorm(1,
c(X[i,]%*%c(c[s,1],as.vector(b[s,1:9]),as.vector(c[s,2:3])))+
elections.simul[df$election[i],1,s]+districts.simul[df$state[i],1,s],
sigma.invalid[s]))

P.simulated[i,2,]<-sapply(1:Multilevelregression[[6]], function(s) rnorm(1,
X[i,]%*%c(c[s,4],as.vector(b[s,10:18]),as.vector(c[s,5:6]))+
 elections.simul[df$election[i],2,s]+districts.simul[df$state[i],2,s],
sigma.absenteeism[s]))

P.expected[i,1,]<-sapply(1:Multilevelregression[[6]], function(s) X[i,]%*%c(c[s,1],as.vector(b[s,1:9]),as.vector(c[s,2:3])))

P.expected[i,2,]<-sapply(1:Multilevelregression[[6]], function(s) X[i,]%*%c(c[s,4],as.vector(b[s,10:18]),as.vector(c[s,5:6])))
}

d1.HLM<-sapply(1:Multilevelregression[[6]],function(s) (P.expected[,1,s]-df$Invalid)^2+(P.expected[,2,s]-df$Absenteeism)^2)
d2.HLM<-sapply(1:Multilevelregression[[6]], function(s)(df$Invalid-P.simulated[,1,s])^2+(df$Absenteeism-P.simulated[,2,s])^2)

######################## TABLE A.4 ####################
									#
matrix(c(sum(apply(d1.HLM,1,mean)), 			#
	sum(apply(d1.compositional,1,mean)),		#
	sum(apply(d2.HLM,1,mean)), 				#
	sum(apply(d2.compositional,1,mean))),		#
	nrow=2,ncol=2, 						#
	dimnames=list(c("HLM", "Compositional"),		#
		c("d1", "d2")))					#
									#
######################################################

#Backup table in .csv
write.csv(matrix(c(sum(apply(d1.HLM,1,mean)), 			
	sum(apply(d1.compositional,1,mean)),		
	sum(apply(d2.HLM,1,mean)), 				
	sum(apply(d2.compositional,1,mean))),		
	nrow=2,ncol=2, 						
	dimnames=list(c("HLM", "Compositional"),		
		c("d1", "d2"))), file="TableA4.csv")

################################ TABLE A.5 ####################################
load("Estimates.Table1.raw")
Coefficients.Absenteeism<-cbind(colMeans(Estimates.Table1$Beta.3), HPDinterval(as.mcmc(Estimates.Table1$Beta.3), p=0.9))
colnames(Coefficients.Absenteeism)<-c("Posterior mean", "90% CI - Low", "90% CI - High")
Coefficients.Invalid<-cbind(colMeans(Estimates.Table1$Beta.2), HPDinterval(as.mcmc(Estimates.Table1$Beta.2), p=0.9))
colnames(Coefficients.Invalid)<-c("Posterior mean", "90% CI - Low", "90% CI - High")


# Table A.5, Column 1
rnd(Coefficients.Absenteeism[-c(1:17),][c(1:6,8,7),],2)

# Table A.5, Column 2
rnd(Coefficients.Invalid[-c(1:17),][c(1:6,8,7),],2)

#Backup table in .csv
write.csv(rnd(Coefficients.Absenteeism[-c(1:17),][c(1:6,8,7),],2), file="TableA5.Col1.csv")
write.csv(rnd(Coefficients.Invalid[-c(1:17),][c(1:6,8,7),],2), file="TableA5.Col2.csv")

################################ FIGURE A.6 #####################################################################################################################################################################################################################################

# Need to compute "marginal effects" first
df.individual<-read.csv("individual.2002.2010.csv")
state<-as.numeric(factor(df.individual$State), levels=unique(df.individual$State))
election<-as.numeric(factor(df.individual$Year), levels=unique(df.individual$Year))
num.districts<-length(unique(state))
num.elections<-length(unique(election))
Y<-cbind(df.individual$Valid, df.individual$Invalid,df.individual$Absenteeism)
y<-c(sapply(1:nrow(Y), function(i) which(Y[i,]==1)))
n<-length(y)          
X<-cbind(1,   df.individual$Income, df.individual$University,
		df.individual$Secondary,
		df.individual$Political.knowledge, df.individual$Urban,
		(df.individual$Candidates-mean(df.individual$Candidates))/(2*sd(df.individual$Candidates)), 
		df.individual$Political.inefficacy, 
		 df.individual$No.representation, 
		df.individual$Dissatisfaction.democracy, 
		(df.individual$Growth-mean(df.individual$Growth))/(2*sd(df.individual$Growth)),	
		(log(df.individual$Inflation)-mean(log(df.individual$Inflation)))/(2*sd(log(df.individual$Inflation))),
		df.individual$TRE, 
		df.individual$Illiterates, 
		df.individual$Young,
		df.individual$Seniors, 
		(df.individual$Competitiveness-mean(df.individual$Competitiveness))/(2*sd(df.individual$Competitiveness)),
		(df.individual$Age-mean(df.individual$Age))/(2*sd(df.individual$Age)),
		 df.individual$Female, df.individual$Married,
		df.individual$Religion, df.individual$Union, 
		df.individual$PT, df.individual$PSDB, df.individual$PMDB
		)
colnames(X)<-c("Intercept", "Income", "University", "Secondary",
		"Pol.Knowledge",
			"Urban", "Candidates", "Political Inefficacy", "No Political Representation",
		"Dissatisfaction w. Democracy", "Growth", "Inflation",
		"Clearance Rate", "Illiterates", "Young", "Seniors", 
		"Competitiveness", 
	"Age", "Female", "Married", "Religion", "Union",
	 "PT", "PSDB", "PMDB")

# Compute marginal effects
dummies<-c(19:22)
continuous<-c(18)
dummies.partisan<-c(23:25)

Marginal.prob.controls<-matrix(NA, nrow=(length(dummies)+length(continuous)+length(dummies.partisan))*2, ncol=5)
colnames(Marginal.prob.controls)<-c("Dep.Variable", "Covariate", "Predictive Comparison", "90% HPD -Lower", "90% HPD - Higher")
Marginal.prob.controls[,1]<-c(rep("Absenteeism",(length(dummies)+length(dummies.partisan)+length(continuous))) ,rep("Invalid Voting",(length(dummies)+length(dummies.partisan)+length(continuous))))
Marginal.prob.controls[,2]<-rep(colnames(X)[sort(c(dummies, dummies.partisan, continuous))],2)

for (j in dummies) {
X.0<-X
X.0[,j]<-0
X.1<-X
X.1[,j]<-1

Absenteeism<-sapply(1:dim(Estimates.Table1$Beta.2)[1], function(s) mean(
exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])/(1+
exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])) -
exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])/(1+
exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])))
)

Invalid<-sapply(1:dim(Estimates.Table1$Beta.3)[1], function(s) mean(
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])/(1+
exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])) -
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])/(1+
exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])))
)

Marginal.prob.controls[grep(colnames(X)[j],Marginal.prob.controls[,2]),3]<-c(mean(Absenteeism), mean(Invalid))
Marginal.prob.controls[grep(colnames(X)[j],Marginal.prob.controls[,2]),4:5]<-matrix(c(HPDinterval(as.mcmc(Absenteeism), p=0.9), HPDinterval(as.mcmc(Invalid), p=0.9)),ncol=2,byrow=TRUE)
}

for (j in continuous) {
X.0<-X
X.1<-X
X.1[,j]<-X.0[,j]+sd(X[,j])

Absenteeism<-sapply(1:dim(Estimates.Table1$Beta.2)[1], function(s) mean(
exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])/(1+
exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])) -
exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])/(1+
exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])))
)

Invalid<-sapply(1:dim(Estimates.Table1$Beta.3)[1], function(s) mean(
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])/(1+
exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])) -
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])/(1+
exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])))
)

Marginal.prob.controls[grep(colnames(X)[j],Marginal.prob.controls[,2]),3]<-c(mean(Absenteeism), mean(Invalid))
Marginal.prob.controls[grep(colnames(X)[j],Marginal.prob.controls[,2]),4:5]<-matrix(c(HPDinterval(as.mcmc(Absenteeism), p=0.9), HPDinterval(as.mcmc(Invalid), p=0.9)),ncol=2,byrow=TRUE)
}

for (j in dummies.partisan) {
X.0<-X
X.1<-X
X.0[,dummies.partisan]<-0
X.1[,dummies.partisan]<-0
X.1[,j]<-1
X.0[,j]<-0

Absenteeism<-sapply(1:dim(Estimates.Table1$Beta.2)[1], function(s) mean(
exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])/(1+
exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])) -
exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])/(1+
exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])))
)

Invalid<-sapply(1:dim(Estimates.Table1$Beta.3)[1], function(s) mean(
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])/(1+
exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])) -
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])/(1+
exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])+
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1])))
)

Marginal.prob.controls[grep(colnames(X)[j],Marginal.prob.controls[,2]),3]<-c(mean(Absenteeism), mean(Invalid))
Marginal.prob.controls[grep(colnames(X)[j],Marginal.prob.controls[,2]),4:5]<-matrix(c(HPDinterval(as.mcmc(Absenteeism), p=0.9), HPDinterval(as.mcmc(Invalid), p=0.9)),ncol=2,byrow=TRUE)
}

Marginal.prob.controls<-data.frame(Marginal.prob.controls)
names(Marginal.prob.controls)<-c("Dep.Variable", "Covariate", "Mean","Lo.CI", "Hi.CI")
for (k in 3:ncol(Marginal.prob.controls)) {
Marginal.prob.controls[,k]<-as.numeric(levels(Marginal.prob.controls[,k]))[Marginal.prob.controls[,k]]
}
Marginal.prob.controls[,3:5]<-100*Marginal.prob.controls[,3:5]
Marginal.prob.controls$Covariate<- factor(Marginal.prob.controls$Covariate, 
	levels = rev(c("Age", "Female", "Married",
		"Religion", "Union", "PT", "PMDB", "PSDB")),
	labels=rev(c("Age", "Female", "Married",
		"Religion", "Union", "PT", "PMDB", "PSDB")))


p.A6<-ggplot(Marginal.prob.controls,
 aes(y=Covariate, x=Mean)) + 
    geom_errorbarh(aes(xmin=Lo.CI, 
	xmax=Hi.CI, height=.2), colour="black") + facet_wrap(~Dep.Variable)+
    geom_point(size=3)+
theme_bw()+xlab("Marginal effect 
(in percentage points)")+ylab("")+
geom_hline(yintercept=0, linetype="longdash", colour="gray40")+
theme(strip.text.x = element_text(size = 12))+
scale_x_continuous(breaks=c(-20, -10, 0, 10, 20))+
geom_vline(xintercept=0, linetype="longdash", colour="gray40")
p.A6<-p.A6+theme(axis.text.y=element_text(size=12), 
	axis.text.x=element_text(size=11),
	axis.title.x=element_text(vjust=-0.5))+
	theme(panel.grid.major=element_blank(),
	panel.grid.minor=element_blank())+
theme(strip.text.x = element_text(size = 12, face="bold"))

## Figure A.6 - display #
	p.A6

# Save figure A.6
 jpeg("FigureA6.jpg")
 p.A6
 dev.off()

################################ TABLE A.6 #################################################################################################################################################################################################

# Marginal effect on absenteeism versus invalid voting

dummies<-c(2,3,4,9,14,15,16, 19:22)
categorical<-c(5,8,10)
continuous<-c(6,7,11,12,13,17,18)
dummies.partisan<-c(23:25)

Marginal.relative<-matrix(NA, nrow=(length(dummies)+length(continuous)+length(categorical)+length(dummies.partisan)), ncol=4)
colnames(Marginal.relative)<-c("Covariate", "Predictive Comparison", "90% HPD -Lower", "90% HPD - Higher")
Marginal.relative[,1]<-colnames(X)[sort(c(dummies, categorical, dummies.partisan, continuous))]

for (j in dummies) {
X.0<-X
X.0[,j]<-0
X.1<-X
X.1[,j]<-1

Relative<-sapply(1:dim(Estimates.Table1$Beta.2)[1], function(s) mean(
(exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])/
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1]))-
(exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])/
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1]))
)
)

Marginal.relative[grep(colnames(X)[j],Marginal.relative[,1]),2:4]<-c(mean(Relative), quantile(Relative,c(0.05,0.95)))
}


for (j in continuous) {
X.0<-X
X.1<-X
X.1[,j]<-X.0[,j]+sd(X[,j])

Relative<-sapply(1:dim(Estimates.Table1$Beta.2)[1], function(s) mean(
(exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])/
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1]))-
(exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])/
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1]))
)
)

Marginal.relative[grep(colnames(X)[j],Marginal.relative[,1]),2:4]<-c(mean(Relative), quantile(Relative,c(0.05,0.95)))
}

for (j in categorical) {
X.0<-X
X.1<-X
X.1[,j]<-X.0[,j]+1

Relative<-sapply(1:dim(Estimates.Table1$Beta.2)[1], function(s) mean(
(exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])/
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1]))-
(exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])/
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1]))
)
)

Marginal.relative[grep(colnames(X)[j],Marginal.relative[,1]),2:4]<-c(mean(Relative), quantile(Relative,c(0.05,0.95)))
}

for (j in dummies.partisan) {
X.0<-X
X.1<-X
X.0[,dummies.partisan]<-0
X.1[,dummies.partisan]<-0
X.1[,j]<-1
X.0[,j]<-0

Relative<-sapply(1:dim(Estimates.Table1$Beta.2)[1], function(s) mean(
(exp(X.1%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])/
exp(X.1%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1]))-
(exp(X.0%*%Estimates.Table1$Beta.3[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,2]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,2])/
exp(X.0%*%Estimates.Table1$Beta.2[s,]+matrix(Estimates.Table1$Reffect.Election[s,],ncol=2)[election,1]+matrix(Estimates.Table1$Reffect.District[s,],ncol=2)[state,1]))
)
)

Marginal.relative[grep(colnames(X)[j],Marginal.relative[,1]),2:4]<-c(mean(Relative), quantile(Relative,c(0.05,0.95)))
}


Marginal.relative.df<-data.frame(Marginal.relative)
names(Marginal.relative.df)<-c("Covariate", "Mean Difference", "90% CI, Low",
		"90% CI, High")
for (k in 2:ncol(Marginal.relative.df)) {
Marginal.relative.df[,k]<-rnd(as.numeric(levels(Marginal.relative.df[,k]))[Marginal.relative.df[,k]],2)
}
Marginal.relative.df$Covariate<- factor(Marginal.relative.df$Covariate, 
	levels = c("Income", "University", "Secondary",
		"Pol.Knowledge", "Urban", "Candidates",
		"Political Inefficacy", "No Political Representation",
		"Dissatisfaction w. Democracy", "Growth", "Inflation",
		"Clearance Rate", "Illiterates", "Young", "Seniors",
		"Competitiveness","Age", "Female", "Married",
		"Religion", "Union", "PT", "PSDB", "PMDB"),
	labels=c("Income", "Education: University", "Education: Secondary",
		 "Political Knowledge", "Urban", "Candidates", 
		 "Political Inefficacy", "No Political Representation",
		 "Dissatisfaction with Democracy", 	"Growth", "Inflation", 
		"Clearance Rate",	"Illiterates", "Young", "Seniors",
		"Competitiveness", "(Log) Age", "Female", "Married",
		"Religion", "Union", "PT", "PSDB", "PMDB"))



########### TABLE A.6 ##############################
     				   				   #
   Marginal.relative.df[c(2,3,1, 4:22, 24,23),]    #
				   				   #
####################################################

#Backup table in .csv
write.csv( Marginal.relative.df[c(2,3,1, 4:22, 24,23),], file="TableA6.csv")

################################### TABLE A.7 #################################################################################################################################################################################################
df.individual<-read.csv("individual.2002.2010.csv")

state<-as.numeric(factor(df.individual$State), levels=unique(df.individual$State))
election<-as.numeric(factor(df.individual$Year), levels=unique(df.individual$Year))
num.districts<-length(unique(state))
num.elections<-length(unique(election))

Y<-cbind(df.individual$Valid, df.individual$Invalid,df.individual$Absenteeism)
y<-c(sapply(1:nrow(Y), function(i) which(Y[i,]==1)))
n<-length(y)          # Number of obs
K<-3			    # Number of outcomes

Matrix.Election<-as.spam(matrix(0,n,num.elections))           
for (i in 1:num.elections) Matrix.Election[election==i,i]<-1
Matrix.District<-as.spam(matrix(0,n,num.districts))
for (i in 1:num.districts) Matrix.District[state==i,i]<-1

X<-cbind(1,   df.individual$Income, df.individual$University,
		df.individual$Secondary,df.individual$Political.knowledge, 
		(df.individual$Candidates-mean(df.individual$Candidates))/(2*sd(df.individual$Candidates)), 
   		(df.individual$Candidates-mean(df.individual$Candidates))/(2*sd(df.individual$Candidates))*df.individual$Income,
   		(df.individual$Candidates-mean(df.individual$Candidates))/(2*sd(df.individual$Candidates))*df.individual$University,
   		(df.individual$Candidates-mean(df.individual$Candidates))/(2*sd(df.individual$Candidates))*df.individual$Secondary,
		(df.individual$Candidates-mean(df.individual$Candidates))/(2*sd(df.individual$Candidates))*df.individual$Political.knowledge,
		df.individual$Political.inefficacy, 
		 df.individual$No.representation, 
		df.individual$Dissatisfaction.democracy, 
		(df.individual$Growth-mean(df.individual$Growth))/(2*sd(df.individual$Growth)),	
		(df.individual$Growth-mean(df.individual$Growth))/(2*sd(df.individual$Growth))*df.individual$Political.inefficacy,
		(df.individual$Growth-mean(df.individual$Growth))/(2*sd(df.individual$Growth))*df.individual$No.representation,
		(df.individual$Growth-mean(df.individual$Growth))/(2*sd(df.individual$Growth))*df.individual$Dissatisfaction.democracy,
		(log(df.individual$Inflation)-mean(log(df.individual$Inflation)))/(2*sd(log(df.individual$Inflation))),
		(log(df.individual$Inflation)-mean(log(df.individual$Inflation)))/(2*sd(log(df.individual$Inflation)))*df.individual$Political.inefficacy ,
		(log(df.individual$Inflation)-mean(log(df.individual$Inflation)))/(2*sd(log(df.individual$Inflation)))*df.individual$No.representation,
		(log(df.individual$Inflation)-mean(log(df.individual$Inflation)))/(2*sd(log(df.individual$Inflation)))*df.individual$Dissatisfaction.democracy,
		 df.individual$Urban, 
		df.individual$TRE, 
		df.individual$Illiterates, 
		df.individual$Young,
		df.individual$Seniors, 
		(df.individual$Competitiveness-mean(df.individual$Competitiveness))/(2*sd(df.individual$Competitiveness)),
		(df.individual$Age-mean(df.individual$Age))/(2*sd(df.individual$Age)),
		 df.individual$Female, df.individual$Married,
		df.individual$Religion, df.individual$Union, 
		df.individual$PT, df.individual$PSDB, df.individual$PMDB
		)
	  
colnames(X)<-c("Intercept", "Income", "University", "Secondary","Pol.Knowledge",
		"Candidates", "Candidates*Income", "Candidates*University", 
		"Candidates*Secondary", "Candidates*Pol.Knowledge",
		"Political Inefficacy", "No Political Representation","Dissatisfaction w. Democracy", 
		"Growth","Growth*Inefficacy" , "Growth*No.Representation", "Growth*Dissatisfaction.Democracy",
		"Inflation","Inflation*Inefficacy" , "Inflation*No.Representation", "Inflation*Dissatisfaction.Democracy",
		"Urban", "Clearance Rate", "Illiterates", "Young", "Seniors", "Competitiveness", 
		"Age", "Female", "Married", "Religion", "Union", "PT", "PSDB", "PMDB")


cost.vars<-2:5
protest.vars<-11:13

# Iterations &  burn-in
 nsim<-2000000				# Number of MCMC Iterations
 thin<-150                          # Thinning Interval (should be a divisor of nsim/2)
 burn<-500000                       # Burn-in
 lastit<-(nsim-burn)/thin	      # Last stored value

# Priors
set.seed(112215)
mu0<-rep(0,ncol(X))
Tau0<-diag(1,ncol(X))
Sigma0<-solve(Tau0)

# Store Values
 Beta2<-Beta3<-array(0,c(lastit,ncol(X),3)) 
 Sigma.election<-array(0, c(lastit,4,3))
 Sigma.district<-array(0, c(lastit,4,3))
 Reffect.district<-array(0,c(lastit,num.districts*2,3))
 Reffect.election<-array(0,c(lastit,num.elections*2,3))
 Sigma.slopes.cost<-array(0, c(lastit,(length(cost.vars)*2)^2,3))
 Sigma.slopes.protest<-array(0, c(lastit,(length(protest.vars)*2)^2,3))
 Reffect.slopes.cost<-array(0,c(lastit,num.districts*2*length(cost.vars),3))
 Reffect.slopes.protest<-array(0,c(lastit,num.elections*2*length(protest.vars),3))

colnames(Beta2)<-colnames(X)
colnames(Beta3)<-colnames(X)

# Start of MCMC algorithm
for (ch in 1:3) {
# Inits
  b2<-rep(0,ncol(X))
  b3<-rep(0,ncol(X))
  cov2<-cov3<-0.01*diag(ncol(X))
  reffect.district<-matrix(0, num.districts,2) 
  sigma.district<-diag(2)
  reffect.election<-matrix(0, num.elections,2) 
  sigma.election<-diag(2)
  reffect.slopes.cost<-matrix(0, num.districts,length(cost.vars)*2)
  reffect.slopes.protest<-matrix(0, num.elections,length(protest.vars)*2)
  sigma.slopes.cost<-diag(2*length(cost.vars))
  sigma.slopes.protest<-diag(2*length(protest.vars))

 A<-0

for (i in 1:nsim) {
 # Calculate Likelihood based on old value of b2 and b3
  eta2<-X%*%b2+reffect.district[state,1]+reffect.election[election,1]+
	rowSums(X[,cost.vars]*reffect.slopes.cost[state,1:length(cost.vars)])+
	rowSums(X[,protest.vars]*reffect.slopes.protest[election,1:length(protest.vars)])
  eta3<-X%*%b3+reffect.district[state,2]+reffect.election[election,2]+
	rowSums(X[,cost.vars]*reffect.slopes.cost[state,(length(cost.vars)+1):(2*length(cost.vars))])+
	rowSums(X[,protest.vars]*reffect.slopes.protest[election,(length(protest.vars)+1):(2*length(protest.vars))])

  p1<-1/(1+exp(eta2)+exp(eta3))
  p2<-p1*exp(eta2)
  p3<-p1*exp(eta3)
  lold<-sum(log(p1)*(y==1)+log(p2)*(y==2)+log(p3)*(y==3))

 # Draw Candidates
  b2new<-b2 + rmvt(1,sigma=.01*cov2,3)  # Tuning parameter to reach appropriate acceptance rate	
  b3new<-b3 + rmvt(1,sigma=.01*cov3,3)   	

 eta2<-X%*%c(b2new)+reffect.district[state,1]+reffect.election[election,1]+
		rowSums(X[,cost.vars]*reffect.slopes.cost[state,1:length(cost.vars)])+
	rowSums(X[,protest.vars]*reffect.slopes.protest[election,1:length(protest.vars)])
  eta3<-X%*%c(b3new)+reffect.district[state,2]+reffect.election[election,2]+
		rowSums(X[,cost.vars]*reffect.slopes.cost[state,(length(cost.vars)+1):(2*length(cost.vars))])+
	rowSums(X[,protest.vars]*reffect.slopes.protest[election,(length(protest.vars)+1):(2*length(protest.vars))])
  p1<-1/(1+exp(eta2)+exp(eta3))
  p2<-p1*exp(eta2)
  p3<-p1*exp(eta3)
  lnew<-sum(log(p1)*(y==1)+log(p2)*(y==2)+log(p3)*(y==3))


# Acceptance prob on log scale
   r<-lnew+dmvnorm(b2new,mu0,Sigma0,log=T)+dmvnorm(b3new,mu0,Sigma0,log=T)-	
	(lold+dmvnorm(b2,mu0,Sigma0,log=T)+dmvnorm(b3,mu0,Sigma0,log=T))
   if(log(runif(1))<r){
    b2<-c(b2new)
    b3<-c(b3new)  
    A<-A+1
    }

#############################################
# Update District random intercepts using MH #
 ############################################
 reffect.district.new<-reffect.district+rmvt(num.districts,sigma=0.001*diag(2),3) 
  eta2.old<-X%*%b2+reffect.district[state,1]+reffect.election[election,1]+
	rowSums(X[,cost.vars]*reffect.slopes.cost[state,1:length(cost.vars)])+
	rowSums(X[,protest.vars]*reffect.slopes.protest[election,1:length(protest.vars)])
  eta3.old<-X%*%b3+reffect.district[state,2]+reffect.election[election,2]+
	rowSums(X[,cost.vars]*reffect.slopes.cost[state,(length(cost.vars)+1):(2*length(cost.vars))])+
	rowSums(X[,protest.vars]*reffect.slopes.protest[election,(length(protest.vars)+1):(2*length(protest.vars))])
  p1.old<-1/(1+exp(eta2.old)+exp(eta3.old))
  p2.old<-p1.old*exp(eta2.old)
  p3.old<-p1.old*exp(eta3.old)
 l.old<-log(p1.old)*(y==1)+log(p2.old)*(y==2)+log(p3.old)*(y==3)
 eta2.new<-X%*%b2+reffect.district.new[state,1]+reffect.election[election,1]+
	rowSums(X[,cost.vars]*reffect.slopes.cost[state,1:length(cost.vars)])+
	rowSums(X[,protest.vars]*reffect.slopes.protest[election,1:length(protest.vars)])
 eta3.new<-X%*%b3+reffect.district.new[state,2]+reffect.election[election,2]+
	rowSums(X[,cost.vars]*reffect.slopes.cost[state,(length(cost.vars)+1):(2*length(cost.vars))])+
	rowSums(X[,protest.vars]*reffect.slopes.protest[election,(length(protest.vars)+1):(2*length(protest.vars))])
  p1.new<-1/(1+exp(eta2.new)+exp(eta3.new))
  p2.new<-p1.new*exp(eta2.new)
  p3.new<-p1.new*exp(eta3.new)
 l.new<-log(p1.new)*(y==1)+log(p2.new)*(y==2)+log(p3.new)*(y==3)
 ratio<-(t(Matrix.District)%*%(l.new-l.old)) + 
	dmvnorm(reffect.district.new,rep(0,2),sigma.district, log=TRUE)-
	dmvnorm(reffect.district,rep(0,2),sigma.district, log=TRUE)
 u<-(log(runif(num.districts))<  ratio)*1
  reffect.district[!is.na(u) & u==1,]<-reffect.district.new[!is.na(u) & u==1,]
   reffect.district[,1]<-reffect.district[,1]-mean(reffect.district[,1])
   reffect.district[,2]<-reffect.district[,2]-mean(reffect.district[,2])

 ##########################################
 # Update sigma.district random intercept #
 #########################################
sigma.district<-riwish(num.districts+3, .01*diag(2)+crossprod(reffect.district))
    
#############################################
# Update Election random intercepts using MH #
 ############################################
reffect.election.new<-reffect.election+rmvt(num.elections,sigma=0.001*diag(2),3) 
  eta2.old<-X%*%b2+reffect.district[state,1]+reffect.election[election,1]+
	rowSums(X[,cost.vars]*reffect.slopes.cost[state,1:length(cost.vars)])+
	rowSums(X[,protest.vars]*reffect.slopes.protest[election,1:length(protest.vars)])
  eta3.old<-X%*%b3+reffect.district[state,2]+reffect.election[election,2]+
	rowSums(X[,cost.vars]*reffect.slopes.cost[state,(length(cost.vars)+1):(2*length(cost.vars))])+
	rowSums(X[,protest.vars]*reffect.slopes.protest[election,(length(protest.vars)+1):(2*length(protest.vars))])
  p1.old<-1/(1+exp(eta2.old)+exp(eta3.old))
  p2.old<-p1.old*exp(eta2.old)
  p3.old<-p1.old*exp(eta3.old)
 l.old<-log(p1.old)*(y==1)+log(p2.old)*(y==2)+log(p3.old)*(y==3)
 eta2.new<-X%*%b2+reffect.district[state,1]+reffect.election.new[election,1]+
	rowSums(X[,cost.vars]*reffect.slopes.cost[state,1:length(cost.vars)])+
	rowSums(X[,protest.vars]*reffect.slopes.protest[election,1:length(protest.vars)])
 eta3.new<-X%*%b3+reffect.district[state,2]+reffect.election.new[election,2]+
	rowSums(X[,cost.vars]*reffect.slopes.cost[state,(length(cost.vars)+1):(2*length(cost.vars))])+
	rowSums(X[,protest.vars]*reffect.slopes.protest[election,(length(protest.vars)+1):(2*length(protest.vars))])
  p1.new<-1/(1+exp(eta2.new)+exp(eta3.new))
  p2.new<-p1.new*exp(eta2.new)
  p3.new<-p1.new*exp(eta3.new)
 l.new<-log(p1.new)*(y==1)+log(p2.new)*(y==2)+log(p3.new)*(y==3)
 ratio<-(t(Matrix.Election)%*%(l.new-l.old)) + 
	dmvnorm(reffect.election.new,rep(0,2),sigma.election, log=TRUE)-
	dmvnorm(reffect.election,rep(0,2),sigma.election, log=TRUE)
 u<-(log(runif(num.elections))<  ratio)*1
  reffect.election[!is.na(u) & u==1,]<-reffect.election.new[!is.na(u) & u==1,]
   reffect.election[,1]<-reffect.election[,1]-mean(reffect.election[,1])
   reffect.election[,2]<-reffect.election[,2]-mean(reffect.election[,2])

 #########################################
 # Update sigma.election random intercepts#
 #########################################
sigma.election<-riwish(num.elections+3, diag(2)+crossprod(reffect.election))

##################################################
# Update random slopes of cost variables using MH #
 ##################################################
reffect.slopes.cost.new<-reffect.slopes.cost+rmvt(num.districts,sigma=0.001*diag(2*length(cost.vars)),3) 
 eta2.old<-X%*%b2+reffect.district[state,1]+reffect.election[election,1]+
	rowSums(X[,cost.vars]*reffect.slopes.cost[state,1:length(cost.vars)])+
	rowSums(X[,protest.vars]*reffect.slopes.protest[election,1:length(protest.vars)])
  eta3.old<-X%*%b3+reffect.district[state,2]+reffect.election[election,2]+
	rowSums(X[,cost.vars]*reffect.slopes.cost[state,(length(cost.vars)+1):(2*length(cost.vars))])+
	rowSums(X[,protest.vars]*reffect.slopes.protest[election,(length(protest.vars)+1):(2*length(protest.vars))])
  p1.old<-1/(1+exp(eta2.old)+exp(eta3.old))
  p2.old<-p1.old*exp(eta2.old)
  p3.old<-p1.old*exp(eta3.old)
 l.old<-log(p1.old)*(y==1)+log(p2.old)*(y==2)+log(p3.old)*(y==3)
 eta2.new<-X%*%b2+reffect.district[state,1]+reffect.election[election,1]+
	rowSums(X[,cost.vars]*reffect.slopes.cost.new[state,1:length(cost.vars)])+
	rowSums(X[,protest.vars]*reffect.slopes.protest[election,1:length(protest.vars)])
 eta3.new<-X%*%b3+reffect.district[state,2]+reffect.election[election,2]+
	rowSums(X[,cost.vars]*reffect.slopes.cost.new[state,(length(cost.vars)+1):(2*length(cost.vars))])+
	rowSums(X[,protest.vars]*reffect.slopes.protest[election,(length(protest.vars)+1):(2*length(protest.vars))])
  p1.new<-1/(1+exp(eta2.new)+exp(eta3.new))
  p2.new<-p1.new*exp(eta2.new)
  p3.new<-p1.new*exp(eta3.new)
 l.new<-log(p1.new)*(y==1)+log(p2.new)*(y==2)+log(p3.new)*(y==3)
 ratio<-(t(Matrix.District)%*%(l.new-l.old)) + 
	dmvnorm(reffect.slopes.cost.new,rep(0,2*length(cost.vars)),sigma.slopes.cost, log=TRUE)-
	dmvnorm(reffect.slopes.cost,rep(0,2*length(cost.vars)),sigma.slopes.cost, log=TRUE)
 u<-(log(runif(num.districts))< ratio)*1
  reffect.slopes.cost[!is.na(u) & u==1,]<-reffect.slopes.cost.new[!is.na(u) & u==1,]
   reffect.slopes.cost[,1]<-reffect.slopes.cost[,1]-mean(reffect.slopes.cost[,1])
   reffect.slopes.cost[,2]<-reffect.slopes.cost[,2]-mean(reffect.slopes.cost[,2])

 ##########################################
 # Update sigma.slopes.cost			 #
 #########################################
sigma.slopes.cost<-riwish(num.districts+(2*length(cost.vars)+1), .01*diag(2*length(cost.vars))+crossprod(reffect.slopes.cost))
    
##################################################
# Update random slopes of protest variables using MH #
 ##################################################
reffect.slopes.protest.new<-reffect.slopes.protest+rmvt(num.elections,sigma=0.001*diag(2*length(protest.vars)),3) 
  eta2.old<-X%*%b2+reffect.district[state,1]+reffect.election[election,1]+
	rowSums(X[,cost.vars]*reffect.slopes.cost[state,1:length(cost.vars)])+
	rowSums(X[,protest.vars]*reffect.slopes.protest[election,1:length(protest.vars)])
  eta3.old<-X%*%b3+reffect.district[state,2]+reffect.election[election,2]+
	rowSums(X[,cost.vars]*reffect.slopes.cost[state,(length(cost.vars)+1):(2*length(cost.vars))])+
	rowSums(X[,protest.vars]*reffect.slopes.protest[election,(length(protest.vars)+1):(2*length(protest.vars))])
  p1.old<-1/(1+exp(eta2.old)+exp(eta3.old))
  p2.old<-p1.old*exp(eta2.old)
  p3.old<-p1.old*exp(eta3.old)
 l.old<-log(p1.old)*(y==1)+log(p2.old)*(y==2)+log(p3.old)*(y==3)
 eta2.new<-X%*%b2+reffect.district[state,1]+reffect.election[election,1]+
	rowSums(X[,cost.vars]*reffect.slopes.cost[state,1:length(cost.vars)])+
	rowSums(X[,protest.vars]*reffect.slopes.protest.new[election,1:length(protest.vars)])
 eta3.new<-X%*%b3+reffect.district[state,2]+reffect.election[election,2]+
	rowSums(X[,cost.vars]*reffect.slopes.cost[state,(length(cost.vars)+1):(2*length(cost.vars))])+
	rowSums(X[,protest.vars]*reffect.slopes.protest.new[election,(length(protest.vars)+1):(2*length(protest.vars))])
  p1.new<-1/(1+exp(eta2.new)+exp(eta3.new))
  p2.new<-p1.new*exp(eta2.new)
  p3.new<-p1.new*exp(eta3.new)
 l.new<-log(p1.new)*(y==1)+log(p2.new)*(y==2)+log(p3.new)*(y==3)
 ratio<-(t(Matrix.Election)%*%(l.new-l.old)) + 
	dmvnorm(reffect.slopes.protest.new,rep(0,2*length(protest.vars)),sigma.slopes.protest, log=TRUE)-
	dmvnorm(reffect.slopes.protest,rep(0,2*length(protest.vars)),sigma.slopes.protest, log=TRUE)
 u<-(log(runif(num.elections))<  ratio)*1
  reffect.slopes.protest[!is.na(u) & u==1,]<-reffect.slopes.protest.new[!is.na(u) & u==1,]
   reffect.slopes.protest[,1]<-reffect.slopes.protest[,1]-mean(reffect.slopes.protest[,1])
   reffect.slopes.protest[,2]<-reffect.slopes.protest[,2]-mean(reffect.slopes.protest[,2])


 ##########################################
 # Update sigma.slopes.protest			 #
 #########################################
sigma.slopes.protest<-riwish(num.elections+(2*length(protest.vars)+1), diag(2*length(protest.vars))+crossprod(reffect.slopes.protest))


# Store Results
  if (i> burn & i%%thin==0) {
    j<-(i-burn)/thin
   Beta2[j,,ch]<-b2
   Beta3[j,,ch]<-b3
   Sigma.district[j,,ch]<-c(sigma.district)
   Sigma.election[j,,ch]<-c(sigma.election)
   Reffect.election[j,,ch]<-as.vector(reffect.election)
   Reffect.district[j,,ch]<-as.vector(reffect.district)
   Sigma.slopes.cost[j,,ch]<-c(sigma.slopes.cost)
   Sigma.slopes.protest[j,,ch]<-c(sigma.slopes.protest)
   Reffect.slopes.cost[j,,ch]<-as.vector(reffect.slopes.cost)
   Reffect.slopes.protest[j,,ch]<-as.vector(reffect.slopes.protest)
}

if (i%%500000==0)  print(c(paste("Chain=",ch), paste("Iteration=",i)))

# Backup saving (if something goes wrong)
if (i==nsim) {

Results<-list(Beta2=Beta2[,,ch], Beta3=Beta3[,,ch], Sigma.district=Sigma.district[,,ch], 
			Sigma.election=Sigma.election[,,ch], 
			Sigma.slopes.cost=Sigma.slopes.cost[,,ch],
			Sigma.slopes.protest=Sigma.slopes.protest[,,ch],
			Reffect.election=Reffect.election[,,ch],
			Reffect.district=Reffect.district[,,ch],
			Reffect.slopes.cost=Reffect.slopes.cost[,,ch],
			Reffect.slopes.protest=Reffect.slopes.protest[,,ch],
			vars=colnames(Beta2))
save(Results, file=paste("Individual-level model with random slopes, chain=", ch))
}

}

} 
# End of MCMC algorithm

# Posterior summaries based on 1,000 convergent samples from each chain
Beta.2<-rbind(Beta2[9001:10000,,1], Beta2[9001:10000,,2], Beta2[9001:10000,,3])
Beta.3<-rbind(Beta3[9001:10000,,1], Beta3[9001:10000,,2], Beta3[9001:10000,,3])
Reffect.Election<-rbind(Reffect.election[9001:10000,,1],Reffect.election[9001:10000,,2],Reffect.election[9001:10000,,3])
Reffect.District<-rbind(Reffect.district[9001:10000,,1],Reffect.district[9001:10000,,2],Reffect.district[9001:10000,,3])
Reffect.Slopes.Cost<-rbind(Reffect.slopes.cost[9001:10000,,1], Reffect.slopes.cost[9001:10000,,2],
			Reffect.slopes.cost[9001:10000,,3])
Reffect.Slopes.Protest<-rbind(Reffect.slopes.protest[9001:10000,,1], Reffect.slopes.protest[9001:10000,,2],
			Reffect.slopes.protest[9001:10000,,3])

Coefficients.Absenteeism<-cbind(colMeans(Beta.3), HPDinterval(as.mcmc(Beta.3), p=0.9))
colnames(Coefficients.Absenteeism)<-c("Posterior mean", "90% CI - Low", "90% CI - High")
Coefficients.Invalid<-cbind(colMeans(Beta.2), HPDinterval(as.mcmc(Beta.2), p=0.9))
colnames(Coefficients.Invalid)<-c("Posterior mean", "90% CI - Low", "90% CI - High")

TableA7.Col1<-rnd(Coefficients.Absenteeism,2)[c(1,3,4,2,5,6,8:9,7,10, 11:33,35,34),]
TableA7.Col2<-rnd(Coefficients.Invalid,2)[c(1,3,4,2,5,6,8:9,7,10, 11:33,35,34),]

# Percent correctly predicted
mean.predicted.probs<-cbind(1/(1+
exp(X%*%apply(Beta.3,2,mean)+matrix(apply(Reffect.Election,2,mean),ncol=2)[election,2]+matrix(apply(Reffect.District,2,mean),ncol=2)[state,2]+
rowSums(X[,cost.vars]*matrix(apply(Reffect.Slopes.Cost,2,mean),ncol=2*length(cost.vars))[state,(length(cost.vars)+1):(2*length(cost.vars))])+
rowSums(X[,protest.vars]*matrix(apply(Reffect.Slopes.Protest,2,mean),ncol=2*length(protest.vars))[election,(length(protest.vars)+1):(2*length(protest.vars))]))+
exp(X%*%apply(Beta.2,2,mean)+matrix(apply(Reffect.Election,2,mean),ncol=2)[election,1]+matrix(apply(Reffect.District,2,mean),ncol=2)[state,1]+
rowSums(X[,cost.vars]*matrix(apply(Reffect.Slopes.Cost,2,mean),ncol=2*length(cost.vars))[state,1:length(cost.vars)])+
rowSums(X[,protest.vars]*matrix(apply(Reffect.Slopes.Protest,2,mean),ncol=2*length(protest.vars))[election,1:length(protest.vars)]))),
exp(X%*%apply(Beta.2,2,mean)+matrix(apply(Reffect.Election,2,mean),ncol=2)[election,1]+matrix(apply(Reffect.District,2,mean),ncol=2)[state,1]+
rowSums(X[,cost.vars]*matrix(apply(Reffect.Slopes.Cost,2,mean),ncol=2*length(cost.vars))[state,1:length(cost.vars)])+
rowSums(X[,protest.vars]*matrix(apply(Reffect.Slopes.Protest,2,mean),ncol=2*length(protest.vars))[election,1:length(protest.vars)]))/(1+
exp(X%*%apply(Beta.3,2,mean)+matrix(apply(Reffect.Election,2,mean),ncol=2)[election,2]+matrix(apply(Reffect.District,2,mean),ncol=2)[state,2]+
rowSums(X[,cost.vars]*matrix(apply(Reffect.Slopes.Cost,2,mean),ncol=2*length(cost.vars))[state,(length(cost.vars)+1):(2*length(cost.vars))])+
rowSums(X[,protest.vars]*matrix(apply(Reffect.Slopes.Protest,2,mean),ncol=2*length(protest.vars))[election,(length(protest.vars)+1):(2*length(protest.vars))]))+
exp(X%*%apply(Beta.2,2,mean)+matrix(apply(Reffect.Election,2,mean),ncol=2)[election,1]+matrix(apply(Reffect.District,2,mean),ncol=2)[state,1]+
rowSums(X[,cost.vars]*matrix(apply(Reffect.Slopes.Cost,2,mean),ncol=2*length(cost.vars))[state,1:length(cost.vars)])+
rowSums(X[,protest.vars]*matrix(apply(Reffect.Slopes.Protest,2,mean),ncol=2*length(protest.vars))[election,1:length(protest.vars)]))),
exp(X%*%apply(Beta.3,2,mean)+matrix(apply(Reffect.Election,2,mean),ncol=2)[election,2]+matrix(apply(Reffect.District,2,mean),ncol=2)[state,2]+
rowSums(X[,cost.vars]*matrix(apply(Reffect.Slopes.Cost,2,mean),ncol=2*length(cost.vars))[state,(length(cost.vars)+1):(2*length(cost.vars))])+
rowSums(X[,protest.vars]*matrix(apply(Reffect.Slopes.Protest,2,mean),ncol=2*length(protest.vars))[election,(length(protest.vars)+1):(2*length(protest.vars))]))/(1+
exp(X%*%apply(Beta.3,2,mean)+matrix(apply(Reffect.Election,2,mean),ncol=2)[election,2]+matrix(apply(Reffect.District,2,mean),ncol=2)[state,2]+
rowSums(X[,cost.vars]*matrix(apply(Reffect.Slopes.Cost,2,mean),ncol=2*length(cost.vars))[state,(length(cost.vars)+1):(2*length(cost.vars))])+
rowSums(X[,protest.vars]*matrix(apply(Reffect.Slopes.Protest,2,mean),ncol=2*length(protest.vars))[election,(length(protest.vars)+1):(2*length(protest.vars))]))+
exp(X%*%apply(Beta.2,2,mean)+matrix(apply(Reffect.Election,2,mean),ncol=2)[election,1]+matrix(apply(Reffect.District,2,mean),ncol=2)[state,1]+
rowSums(X[,cost.vars]*matrix(apply(Reffect.Slopes.Cost,2,mean),ncol=2*length(cost.vars))[state,1:length(cost.vars)])+
rowSums(X[,protest.vars]*matrix(apply(Reffect.Slopes.Protest,2,mean),ncol=2*length(protest.vars))[election,1:length(protest.vars)])))
)
max.predicted.probs<-apply(mean.predicted.probs,1,max)
b<-sapply(1:nrow(X), function(i) which(mean.predicted.probs[i,]==max.predicted.probs[i]))
c<-sapply(1:nrow(X), function(i) which(Y[i,]==1))

# Chi-squared goodness of fit test 
# H0: Data is consistent with the fitted model
observed<-sapply(1:K, function(k) sum(Y[,k]))
expected<-sapply(1:K, function(k) sum(mean.predicted.probs[,k]))/nrow(X)


########################### ESTIMATES TABLE A.7 #########################
												#
#Column 1											#
TableA7.Col1										#
												#
#Column 2											#
TableA7.Col2										#
												#
# Percent correctly predicted								#
sum(b==c)/length(b)*100									#
												#
#Chi-squared test										#
chisq.test(observed, p=expected)$p.value 							#
												#		#
# Obs												#
nrow(Y)											#
#########################################################################


#Backup table in .csv
write.csv(TableA7.Col1, file="TableA7.Col1.csv")
write.csv(TableA7.Col2, file="TableA7.Col2.csv")


##################################### FIGURE A.7 #################################################################################################################################################################################################
# Need to compute "marginal effects" first
df<-read.csv("district.1986.2014.csv")
log_invalid<-log(df$Invalid/df$Valid)
log_absenteeism<-log(df$Absenteeism/df$Valid)
Y<-cbind(log_invalid, log_absenteeism)
n<-nrow(Y)
df$state<-as.numeric(factor(df$State), levels=unique(df$State))
df$election<-as.numeric(factor(df$Year), levels=unique(df$Year))
num.districts<-length(unique(df$state))
num.elections<-length(unique(df$election))
X<-cbind(1,	df$Urban,df$Income/5,(df$Candidates-mean(df$Candidates))/(2*sd(df$Candidates)),
	df$Illiterates,df$Young,df$Seniors,df$TRE,df$Electronic.Vote,
	(df$Competitiveness-mean(df$Competitiveness))/(2*sd(df$Competitiveness)),
	(df$Growth-mean(df$Growth))/(2*sd(df$Growth)),	
	(log(df$Inflation)-mean(log(df$Inflation)))/(2*sd(log(df$Inflation))))
colnames(X)<-c("Intercept", "Urban", "Income", "Candidates","Illiterates", "Young",
			"Seniors", "Clearance Rate", "Electronic Vote",
			"Competitiveness",
			"Growth", "Inflation")
n<-nrow(X)

load("Compositional.Cols12.raw")
attach.bugs(Compositional.Cols12, overwrite=TRUE)

covars<-c(2,4:9)
Marginal.FA7<-matrix(NA, nrow=length(covars)*2, ncol=5)
colnames(Marginal.FA7)<-c("Dep.Variable", "Covariate", "Mean",
				"90% CI - Lower", "90% CI - Higher")
Marginal.FA7[,1]<-c(rep("Absenteeism", length(covars)),rep("Invalid Voting",length(covars)))
Marginal.FA7[,2]<-rep(colnames(X)[covars],2)

for (j in covars) {
X.effect<-X
X.effect[,j]<-X[,j]+sd(X[,j])

Effect.Absenteeism<-sapply(1:Compositional.Cols12[[6]],function(s)
	mean(exp(X.effect%*%c(c[s,4],b[s,10:18],c[s,5:6])+
	 alpha[s,df$election,2]+beta[s,df$state,2])-
	exp(X%*%c(c[s,4],b[s,10:18],c[s,5:6])+
	alpha[s,df$election,2]+beta[s,df$state,2])))

Effect.Invalid<-sapply(1:Compositional.Cols12[[6]],function(s)
	mean(exp(X.effect%*%c(c[s,1],b[s,1:9],c[s,2:3])+
	 alpha[s,df$election,1]+beta[s,df$state,1])-
	exp(X%*%c(c[s,1],b[s,1:9],c[s,2:3])+
	alpha[s,df$election,1]+beta[s,df$state,1])))

Marginal.FA7[grep(colnames(X)[j], Marginal.FA7[,2]),3]<-c(mean(Effect.Absenteeism),mean(Effect.Invalid))*100
Marginal.FA7[grep(colnames(X)[j], Marginal.FA7[,2]),4:5]<-matrix(c(HPDinterval(as.mcmc(Effect.Absenteeism), p=0.9)*100, HPDinterval(as.mcmc(Effect.Invalid)*100, p=0.9)),ncol=2,byrow=TRUE)
}

Marginal.FA7<-data.frame(Marginal.FA7)
names(Marginal.FA7)<-c("Dep.Variable", "Covariate", "Mean", "Lo", "Hi")
for (k in 3:ncol(Marginal.FA7)){
Marginal.FA7[,k]<-rnd(as.numeric(levels(Marginal.FA7[,k]))[Marginal.FA7[,k]],2)
}
Marginal.FA7$Covariate<- factor(Marginal.FA7$Covariate, 
	levels = rev(c("Urban", "Candidates",
		"Electronic Vote", "Clearance Rate", "Illiterates",
		"Young", "Seniors")))

p.A7<-ggplot(Marginal.FA7,
 aes(y=Covariate, x=Mean)) + 
    geom_errorbarh(aes(xmin=Lo, 
	xmax=Hi, height=.2), colour="black") + facet_wrap(~Dep.Variable)+
    geom_point(size=3)+
theme_bw()+xlab("Marginal effect 
(in percentage points)")+ylab("")+
geom_hline(yintercept=0, linetype="longdash", colour="gray40")+
theme(strip.text.x = element_text(size = 12))+
scale_x_continuous(breaks=c(-15,-10,-5,0,5))+
geom_vline(xintercept=0, linetype="longdash", colour="gray40")
p.A7<-p.A7+theme(axis.text.y=element_text(size=12), axis.text.x=element_text(size=11),
	axis.title.x=element_text(vjust=-0.5))+
	theme(panel.grid.major=element_blank(),
	panel.grid.minor=element_blank())

## Figure A.7 - display #
	p.A7

# Save figure A.7
 jpeg("FigureA7.jpg")
 p.A7
 dev.off()


################################ TABLE A.8 #################################################################################################################################################################################################
# Estimation of the distrcict-level model for elections from 2002 onwards
df<-read.csv("district.1986.2014.csv")
df.post2002<-subset(df, Year>=2002)

log_invalid<-log(df.post2002$Invalid/df.post2002$Valid)
log_absenteeism<-log(df.post2002$Absenteeism/df.post2002$Valid)
Y<-cbind(log_invalid, log_absenteeism)
n<-nrow(Y)
df.post2002$state<-as.numeric(factor(df.post2002$State), levels=unique(df.post2002$State))
df.post2002$election<-as.numeric(factor(df.post2002$Year), levels=unique(df.post2002$Year))
num.districts<-length(unique(df.post2002$state))
num.elections<-length(unique(df.post2002$election))

k<-2# Number of outcome variables
# Priors 
R<-diag(2)
G0<-diag(2)
H0<-diag(2)
B<-diag(18)
b0=rep(0,18)
C<-diag(6)
c0=rep(0,6)

Compositional.TableA8.data<-list("num.districts"=num.districts,"num.elections"=num.elections,
	"n"=n, "B"=B, "C"=C, "b0"=b0, "c0"=c0,
 	"State"=df.post2002$state, "G0"=G0, "H0"=H0, "R"=R,
	 "Y"=Y, 
     "Urban"=df.post2002$Urban,
	"Income"=df.post2002$Income/5,
	"Candidates"=(df.post2002$Candidates-mean(df.post2002$Candidates))/(2*sd(df.post2002$Candidates)),
	"Illiterate"=df.post2002$Illiterates, 
	 "Young"=df.post2002$Young, 
		"Old"=df.post2002$Seniors, 
	  "TRE"=df.post2002$TRE,
	 "Evote"=df.post2002$Electronic.Vote, 
	 "Competitiveness"=(df.post2002$Competitiveness-mean(df.post2002$Competitiveness))/(2*sd(df.post2002$Competitiveness)),
      "Election"=df.post2002$election,  
	"Growth"=(df.post2002$Growth-mean(df.post2002$Growth))/(2*sd(df.post2002$Growth)),
	"Inflation"=(log(df.post2002$Inflation)-mean(log(df.post2002$Inflation)))/(2*sd(log(df.post2002$Inflation)))	  
)

set.seed(112215)  
Compositional.TableA8.inits<-function() {
list(T=diag(2), G=diag(2), H=diag(2), 
b=rnorm(18), c=rnorm(6),  
beta=matrix(rnorm(num.districts*2),num.districts,2),
alpha=matrix(rnorm(num.elections*2), num.elections,2),
k.e=round(runif(1,1,5)),
w.e=rgamma(n,1,1))} 

Compositional.TableA8.parameters<-c( "b", "c", "alpha", "beta", "w.e", "nu.e", 
	"Sigma", "Sigma.election", "Sigma.district")

Compositional.TableA8<-bugs(Compositional.TableA8.data,
	Compositional.TableA8.inits,Compositional.TableA8.parameters,
	"compositional.19862014.bug", 
n.iter=25000, n.burnin=5000,  
# Set the appropriate path for bugs.directory
bugs.directory="//isad.isadroot.ex.ac.uk/UOE/User/Desktop/WinBUGS14")
attach.bugs(Compositional.TableA8, overwrite=TRUE)

# Check convergence
sum(Compositional.TableA8[[10]][,8]<1.2)/nrow(Compositional.TableA8[[10]])>0.975 #Neelon's criterion

X<-cbind(1,	df.post2002$Urban,df.post2002$Income/5,(df.post2002$Candidates-mean(df.post2002$Candidates))/(2*sd(df.post2002$Candidates)),
	df.post2002$Illiterates,df.post2002$Young,df.post2002$Seniors,df.post2002$TRE,df.post2002$Electronic.Vote,
	(df.post2002$Competitiveness-mean(df.post2002$Competitiveness))/(2*sd(df.post2002$Competitiveness)),
	(df.post2002$Growth-mean(df.post2002$Growth))/(2*sd(df.post2002$Growth)),	
	(log(df.post2002$Inflation)-mean(log(df.post2002$Inflation)))/(2*sd(log(df.post2002$Inflation))))
colnames(X)<-c("Intercept", "Urban", "Income", "Candidates","Illiterates", "Young",
			"Seniors", "Clearance Rate", "Electronic Vote",
			"Competitiveness",
			"Growth", "Inflation")

# Posterior Summaries
Estimates.2002.2014<-matrix(NA, nrow=(ncol(X))*2, ncol=5)
colnames(Estimates.2002.2014)<-c("Dep.Variable", "Covariate", "Posterior Mean",
	"90% CI, Low", "90% CI")
Estimates.2002.2014[,1]<-c(rep("Absenteeism", (ncol(X))),rep("Invalid Voting",(ncol(X))))
Estimates.2002.2014[,2]<-rep(colnames(X),2)

Estimates.2002.2014[Estimates.2002.2014[,1]=="Absenteeism",3:ncol(Estimates.2002.2014)]<-rbind(
	c(mean(c[,4]), HPDinterval(as.mcmc(c[,4]), p=0.9)),
cbind(apply(b[,10:18], 2,mean),HPDinterval(as.mcmc(b[,10:18]), p=0.9)),
			cbind(apply(c[,5:6], 2,mean),HPDinterval(as.mcmc(c[,5:6]), p=0.9)))

Estimates.2002.2014[Estimates.2002.2014[,1]=="Invalid Voting",3:ncol(Estimates.2002.2014)]<-rbind(
	c(mean(c[,1]), HPDinterval(as.mcmc(c[,1]), p=0.9)),
	cbind(apply(b[,1:9], 2,mean),HPDinterval(as.mcmc(b[,1:9]), p=0.9)),
			cbind(apply(c[,2:3], 2,mean),HPDinterval(as.mcmc(c[,2:3]), p=0.9)))

Estimates.2002.2014.df<-data.frame(Estimates.2002.2014)
names(Estimates.2002.2014.df)<-colnames(Estimates.2002.2014)
for (k in 3:ncol(Estimates.2002.2014.df)){
Estimates.2002.2014.df[,k]<-rnd(as.numeric(levels(Estimates.2002.2014.df[,k]))[Estimates.2002.2014.df[,k]],2)
}

TableA8.Col1<-subset(Estimates.2002.2014.df, Dep.Variable=="Absenteeism")
TableA8.Col2<-subset(Estimates.2002.2014.df, Dep.Variable=="Invalid Voting")

########################### ESTIMATES TABLE A.8 ###########################
												  #
# Column 1 (Absenteeism)								  #
TableA8.Col1[c(1,3,2,4,9,8,5,6,7,11,12,10),]					  #						  
												  #
# Column 2 (Invalid Votig)								  #
TableA8.Col2[c(1,3,2,4,9,8,5,6,7,11,12,10),]					  #	  
											        #
#Observations										  #
n												  #
												  #	
###########################################################################

#Backup table in .csv
write.csv(TableA8.Col1[c(1,3,2,4,9,8,5,6,7,11,12,10),], file="TableaA8.Col1.csv")
write.csv(TableA8.Col2[c(1,3,2,4,9,8,5,6,7,11,12,10),], file="TableA8.Col2.csv")

################################ FIGURE A.8 #################################################################################################################################################################################################

# For left panel of the figure: Use estimates  from Table 2, Cols 3-4
load("Compositional.Cols34.raw")
attach.bugs(Compositional.Cols34, overwrite=TRUE)

Illiterates.proportionEvote<-b[,4]
Illiterates.interaction.proportionEvote<-b[,4]+b[,19]

Left.Panel<-cbind("Proportion of electronic votes 
in the district", 
c("Illiterates before Electronic Vote", "Illiterates after Electronic Vote"),
 rbind(c(mean(Illiterates.proportionEvote), HPDinterval(as.mcmc(Illiterates.proportionEvote),p=0.9)),
c(mean(Illiterates.interaction.proportionEvote), HPDinterval(as.mcmc(Illiterates.interaction.proportionEvote),p=0.9))
)
)

# For right panel of the figure: need to re-estimate model
# but using a dummy for elections from 2002 onwards
df<-read.csv("district.1986.2014.csv")
log_invalid<-log(df$Invalid/df$Valid)
log_absenteeism<-log(df$Absenteeism/df$Valid)
Y<-cbind(log_invalid, log_absenteeism)
n<-nrow(Y)
df$state<-as.numeric(factor(df$State), levels=unique(df$State))
df$election<-as.numeric(factor(df$Year), levels=unique(df$Year))
num.districts<-length(unique(df$state))
num.elections<-length(unique(df$election))
k<-2
R<-diag(2)
G0<-diag(2)
H0<-diag(2)
B<-diag(20)
b0=rep(0,20)
C<-diag(6)
c0=rep(0,6)

FigureA8.Rightpanel.data<-list("num.districts"=num.districts,"num.elections"=num.elections,
	"n"=n, "B"=B, "C"=C, "b0"=b0, "c0"=c0,
 	"State"=df$state, "G0"=G0, "H0"=H0, "R"=R,
	 "Y"=Y, 
     "Urban"=df$Urban,
	"Income"=df$Income/5,
	"Candidates"=(df$Candidates-mean(df$Candidates))/(2*sd(df$Candidates)),
	"Illiterate"=df$Illiterates, 
	"Interaction.Illiterates"=df$Illiterates*ifelse(df$Year>=2002,1,0),
	 "Young"=df$Young, 
		"Old"=df$Seniors, 
	  "TRE"=df$TRE,
	 "Evote"=df$Electronic.Vote, 
	 "Competitiveness"=(df$Competitiveness-mean(df$Competitiveness))/(2*sd(df$Competitiveness)),
      "Election"=df$election,  
	"Growth"=(df$Growth-mean(df$Growth))/(2*sd(df$Growth)),
	"Inflation"=(log(df$Inflation)-mean(log(df$Inflation)))/(2*sd(log(df$Inflation)))	  
)
set.seed(112215)  
FigureA8.Rightpanel.inits<-function() {
list(T=diag(2), G=diag(2), H=diag(2), 
b=rnorm(20), c=rnorm(6),  
beta=matrix(rnorm(num.districts*2),num.districts,2),
alpha=matrix(rnorm(num.elections*2), num.elections,2),
k.e=round(runif(1,1,5)),
w.e=rgamma(n,1,1))} 

FigureA8.Rightpanel.parameters<-c( "b", "c", "alpha", "beta", "w.e", "nu.e", 
	"Sigma", "Sigma.election", "Sigma.district")

FigureA8.Rightpanel<-bugs(FigureA8.Rightpanel.data,
	FigureA8.Rightpanel.inits,FigureA8.Rightpanel.parameters,
	"compositional.19862014.interaction.bug", 
n.iter=15000, n.burnin=5000,  
# Set the appropriate path for bugs.directory
bugs.directory="//isad.isadroot.ex.ac.uk/UOE/User/Desktop/WinBUGS14")
attach.bugs(FigureA8.Rightpanel, overwrite=TRUE)

Illiterates.2002<-b[,4]
Illiterates.interaction.2002<-b[,4]+b[,19]

Right.Panel<-cbind("Replacement of paper ballots 
with electronic voting", c("Illiterates before Electronic Vote", "Illiterates after Electronic Vote"),
 rbind(c(mean(Illiterates.2002), 
HPDinterval(as.mcmc(Illiterates.2002),p=0.9)
),
c(mean(Illiterates.interaction.2002), 
HPDinterval(as.mcmc(Illiterates.interaction.2002),p=0.9)))
)

Interactions<-rbind(Left.Panel, Right.Panel)
Interactions<-data.frame(Interactions)
names(Interactions)<-c("Year", "Coefficient", "Mean", "Lo", "Hi")
for (k in 3:5) {
Interactions[,k]<-as.numeric(levels(Interactions[,k]))[Interactions[,k]] 
}
Interactions$Coefficient<- factor(Interactions$Coefficient, levels = c("Illiterates before Electronic Vote", "Illiterates after Electronic Vote"))

p.A8<-ggplot(Interactions,
 aes(x=Coefficient, y=Mean)) + 
    geom_errorbar(aes(ymin=Lo, 
	ymax=Hi), colour="black", 
	 width=.3) + facet_wrap(~Year)+
    geom_point(size=3)+
scale_x_discrete( breaks=c("Illiterates before Electronic Vote", "Illiterates after Electronic Vote"),
            labels=c("Illiterates before 
		Electronic Vote", "Illiterates after 
		Electronic Vote"))+
scale_y_continuous(breaks=c(-1,-.5,0,0.5,1,1.5,2))+
theme_bw()+xlab("")+ylab("Estimate")+
geom_hline(yintercept=0, linetype="longdash", colour="gray40")+
theme(strip.text.x = element_text(size = 12))+
theme(panel.grid.major=element_blank(),
	panel.grid.minor=element_blank())


## Figure A.8 - display #
	p.A8

# Save figure A.8
 jpeg("FigureA8.jpg")
 p.A8
 dev.off()

################################ TABLE A.9 #################################################################################################################################################################################################
# Need to estimate a new model for this table

df.TableA9<-read.csv("supplementary.district.1974.2014.csv")
log_invalid<-log(df.TableA9$Invalid/df.TableA9$Valid)
log_absenteeism<-log(df.TableA9$Absenteeism/df.TableA9$Valid)
Y<-cbind(log_invalid, log_absenteeism)
n<-nrow(Y)
df.TableA9$state<-as.numeric(factor(df.TableA9$State), levels=unique(df.TableA9$State))
df.TableA9$election<-as.numeric(factor(df.TableA9$Year), levels=unique(df.TableA9$Year))
num.districts<-length(unique(df.TableA9$state))
num.elections<-length(unique(df.TableA9$election))

k<-2
R<-diag(2)
G0<-diag(2)
H0<-diag(2)
B<-diag(16)
b0=rep(0,16)
C<-diag(6)
c0=rep(0,6)

Compositional.TableA9.data<-list("num.districts"=num.districts,"num.elections"=num.elections,
	"n"=n, "B"=B, "C"=C, "b0"=b0, "c0"=c0,
 	"State"=df.TableA9$state, "G0"=G0, "H0"=H0, "R"=R,
	 "Y"=Y, 
     "Urban"=df.TableA9$Urban,
	"Candidates"=(df.TableA9$Candidates-mean(df.TableA9$Candidates))/(2*sd(df.TableA9$Candidates)),
	"Illiterate"=df.TableA9$Illiterates, 
	 "Young"=df.TableA9$Young, 
		"Old"=df.TableA9$Seniors, 
	 "Evote"=df.TableA9$Electronic.Vote, 
      "Election"=df.TableA9$election, 
	 "Competitiveness"=(df.TableA9$Competitiveness-mean(df.TableA9$Competitiveness))/(2*sd(df.TableA9$Competitiveness)),
	"Growth"=(df.TableA9$Growth-mean(df.TableA9$Growth))/(2*sd(df.TableA9$Growth)),
	"Inflation"=(log(df.TableA9$Inflation)-mean(log(df.TableA9$Inflation)))/(2*sd(log(df.TableA9$Inflation))),
	"Free.Elections"=(df.TableA9$Dif.Inverse.Fscores-mean(df.TableA9$Dif.Inverse.Fscores))/(2*sd(df.TableA9$Dif.Inverse.Fscores)) 
)

set.seed(112215)  
Compositional.TableA9.inits<-function() {
list(T=diag(2), G=diag(2), H=diag(2), 
b=rnorm(16), c=rnorm(6),  
beta=matrix(rnorm(num.districts*2),num.districts,2),
alpha=matrix(rnorm(num.elections*2), num.elections,2),
k.e=round(runif(1,1,5)),
w.e=rgamma(n,1,1))} 

Compositional.TableA9.parameters<-c( "b", "c", "alpha", "beta", "w.e", "nu.e", 
	"Sigma", "Sigma.election", "Sigma.district")

Compositional.TableA9<-bugs(Compositional.TableA9.data,
	Compositional.TableA9.inits,Compositional.TableA9.parameters,
	"compositional.longdistrict.bug", 
n.iter=15000, n.burnin=5000, 
# Set the appropriate path for bugs.directory
bugs.directory="//isad.isadroot.ex.ac.uk/UOE/User/Desktop/WinBUGS14")
attach.bugs(Compositional.TableA9, overwrite=TRUE)

# Check convergence
sum(Compositional.TableA9[[10]][,8]<1.2)/nrow(Compositional.TableA9[[10]])>0.975 #Neelon's criterion

X<-cbind(1,	df.TableA9$Urban,
	(df.TableA9$Candidates-mean(df.TableA9$Candidates))/(2*sd(df.TableA9$Candidates)),
	df.TableA9$Illiterates,df.TableA9$Young,df.TableA9$Seniors,df.TableA9$Electronic.Vote, 
	(df.TableA9$Competitiveness-mean(df.TableA9$Competitiveness))/(2*sd(df.TableA9$Competitiveness)),
	(df.TableA9$Dif.Inverse.Fscores-mean(df.TableA9$Dif.Inverse.Fscores))/(2*sd(df.TableA9$Dif.Inverse.Fscores)),
	(df.TableA9$Growth-mean(df.TableA9$Growth))/(2*sd(df.TableA9$Growth)),	
	(log(df.TableA9$Inflation)-mean(log(df.TableA9$Inflation)))/(2*sd(log(df.TableA9$Inflation))))
colnames(X)<-c("Intercept", "Urban", "Candidates","Illiterates", "Young",
			"Seniors","Electronic Vote",
			"Competitiveness",
			"P.Rights & C.Liberties",
			"Growth", "Inflation")

# Posterior Summaries
Estimates.TableA9<-matrix(NA, nrow=(ncol(X)*2), ncol=5)
colnames(Estimates.TableA9)<-c("Dep.Variable", "Covariate", "Posterior Mean",
	"90% CI, Low", "90% CI")
Estimates.TableA9[,1]<-c(rep("Absenteeism", ncol(X)),rep("Invalid Voting",ncol(X)))
Estimates.TableA9[,2]<-rep(colnames(X),2)

Estimates.TableA9[Estimates.TableA9[,1]=="Absenteeism",3:ncol(Estimates.TableA9)]<-rbind(c(mean(c[,4]), HPDinterval(as.mcmc(c[,4]),p=0.9)),
			cbind(apply(b[,9:16], 2,mean),HPDinterval(as.mcmc(b[,9:16]), p=0.9)),
			cbind(apply(c[,5:6], 2,mean),HPDinterval(as.mcmc(c[,5:6]), p=0.9)))

Estimates.TableA9[Estimates.TableA9[,1]=="Invalid Voting",3:ncol(Estimates.TableA9)]<-rbind(c(mean(c[,1]), HPDinterval(as.mcmc(c[,1]),p=0.9)),
		cbind(apply(b[,1:8], 2,mean),HPDinterval(as.mcmc(b[,1:8]), p=0.9)),
			cbind(apply(c[,2:3], 2,mean),HPDinterval(as.mcmc(c[,2:3]), p=0.9)))

Estimates.TableA9.df<-data.frame(Estimates.TableA9)
names(Estimates.TableA9.df)<-colnames(Estimates.TableA9)
for (k in 3:ncol(Estimates.TableA9.df)){
Estimates.TableA9.df[,k]<-rnd(as.numeric(levels(Estimates.TableA9.df[,k]))[Estimates.TableA9.df[,k]],2)
}

TableA9.Col1<-subset(Estimates.TableA9.df, Dep.Variable=="Absenteeism")
TableA9.Col2<-subset(Estimates.TableA9.df, Dep.Variable=="Invalid Voting")

########################### ESTIMATES TABLE A.9 ###########################
												  #
# Column 1 (Absenteeism)								  #				  #
TableA9.Col1[c(1,2,3,7,9,10,11,4,5,6,8),]						  #
												  #
# Column 2 (Invalid Voting)								  #				  #
TableA9.Col2[c(1,2,3,7,9,10,11,4,5,6,8),]						  #
												  #
# Deviance Information Criterion 							  #
rnd(Compositional.TableA9[[10]][dim(Compositional.TableA9[[10]])][1],2)   #
										        	  #
#Observations										  #
n												  #
												  #	
###########################################################################

#Backup table in .csv
write.csv(TableA9.Col1[c(1,2,3,7,9,10,11,4,5,6,8),], file="TableA9.Col1.csv")
write.csv(TableA9.Col2[c(1,2,3,7,9,10,11,4,5,6,8),], file="TableA9.Col2.csv")


################################ FIGURE A.9 #################################################################################################################################################################################################

# First need to compute the impact of each set of variables
# on the probabilities of Absenteeism & Invalid Votig
df<-read.csv("district.1986.2014.csv")
log_invalid<-log(df$Invalid/df$Valid)
log_absenteeism<-log(df$Absenteeism/df$Valid)
Y<-cbind(log_invalid, log_absenteeism)
n<-nrow(Y)
df$state<-as.numeric(factor(df$State), levels=unique(df$State))
df$election<-as.numeric(factor(df$Year), levels=unique(df$Year))
num.districts<-length(unique(df$state))
num.elections<-length(unique(df$election))
X<-cbind(1,	df$Urban,df$Income/5,(df$Candidates-mean(df$Candidates))/(2*sd(df$Candidates)),
	df$Illiterates,df$Young,df$Seniors,df$TRE,df$Electronic.Vote,
	(df$Competitiveness-mean(df$Competitiveness))/(2*sd(df$Competitiveness)),
	(df$Growth-mean(df$Growth))/(2*sd(df$Growth)),	
	(log(df$Inflation)-mean(log(df$Inflation)))/(2*sd(log(df$Inflation))))
colnames(X)<-c("Intercept", "Urban", "Income", "Candidates","Illiterates", "Young",
			"Seniors", "Clearance Rate", "Electronic Vote",
			"Competitiveness",
			"Growth", "Inflation")

load("Compositional.Cols12.raw")
attach.bugs(Compositional.Cols12, overwrite=TRUE)

# Impact of different sets of variables on probabilities of non-voting
X.0<-c(1, max(X[,2]),mean(X[,3]), min(X[,4]),
		min(X[,5]), mean(X[,6]), min(X[,7]),
		max(X[,8]), max(X[,9]), 
		mean(X[,10]), mean(X[,11]), mean(X[,12]))

X.1<-c(1,min(X[,2]), mean(X[,3]), max(X[,4]), 
		max(X[,5]), mean(X[,6]), max(X[,7]),
		min(X[,8]), min(X[,9]),
		mean(X[,10]), mean(X[,11]), mean(X[,12]))
	
mean.alpha.1<-apply(alpha[,,1],1,mean)
mean.alpha.2<-apply(alpha[,,2],1,mean)
mean.beta.1<-apply(beta[,,1],1,mean)
mean.beta.2<-apply(beta[,,2],1,mean)

Absenteeism.aggregate.0<-mean(sapply(1:Compositional.Cols12[[6]],function(s)
mean(
exp(X.0%*%c(c[s,4],b[s,10:18],c[s,5:6])+
	 mean.alpha.2[s]+mean.beta.2[s])/(1+
exp(X.0%*%c(c[s,4],b[s,10:18],c[s,5:6])+
	 mean.alpha.2[s]+mean.beta.2[s])+
exp(X.0%*%c(c[s,1],b[s,1:9],c[s,2:3])+
	 mean.alpha.1[s]+mean.beta.1[s])))
)
)

Absenteeism.aggregate.1<-mean(sapply(1:Compositional.Cols12[[6]],function(s)
mean(
exp(X.1%*%c(c[s,4],b[s,10:18],c[s,5:6])+
	 mean.alpha.2[s]+mean.beta.2[s])/(1+
exp(X.1%*%c(c[s,4],b[s,10:18],c[s,5:6])+
	 mean.alpha.2[s]+mean.beta.2[s])+
exp(X.1%*%c(c[s,1],b[s,1:9],c[s,2:3])+
	 mean.alpha.1[s]+mean.beta.1[s])))
)
)

Invalid.aggregate.0<-mean(sapply(1:Compositional.Cols12[[6]],function(s)
mean(
exp(X.0%*%matrix(c(c[s,1],b[s,1:9],c[s,2:3]))+
	mean.alpha.1[s]+mean.beta.1[s])/(1+
exp(X.0%*%matrix(c(c[s,4],b[s,10:18],c[s,5:6]))+
	 mean.alpha.2[s]+mean.beta.2[s])+
exp(X.0%*%matrix(c(c[s,1],b[s,1:9],c[s,2:3]))+
	mean.alpha.1[s]+mean.beta.1[s])))
)
)

Invalid.aggregate.1<-mean(sapply(1:Compositional.Cols12[[6]],function(s)
mean(
exp(X.1%*%matrix(c(c[s,1],b[s,1:9],c[s,2:3]))+
	mean.alpha.1[s]+mean.beta.1[s])/(1+
exp(X.1%*%matrix(c(c[s,4],b[s,10:18],c[s,5:6]))+
	 mean.alpha.2[s]+mean.beta.2[s])+
exp(X.1%*%matrix(c(c[s,1],b[s,1:9],c[s,2:3]))+
	mean.alpha.1[s]+mean.beta.1[s]))))
)

dd<-rbind(cbind("Absenteeism",c("Minimal-abstention", "Impact of socio-demographic 
	& politico-institutional variables"), 
c(mean(Absenteeism.aggregate.0),
mean(Absenteeism.aggregate.1))),
cbind("Invalid Voting",c("Minimal-abstention", "Impact of socio-demographic 
	& politico-institutional variables"),
c(mean(Invalid.aggregate.0),
mean(Invalid.aggregate.1)))
)
dd<-data.frame(dd)
names(dd)<-c("Dep.Variable", "Scenario", "Probability")
for (k in 3:3) {
dd[,k]<-as.numeric(levels(dd[,k]))[dd[,k]]
}
dd$Scenario<- factor(dd$Scenario, 
	levels = rev(c("Minimal-abstention",
	 "Impact of socio-demographic \n& politico-institutional variables")))

p.A9<-ggplot(dd, aes(x=Scenario, y=Probability))+
	geom_bar(stat="identity", fill="gray60", width=0.1)+
	facet_grid(.~Dep.Variable)+theme_bw()+
	ylab("Proportions of non-voters in the district")+xlab("")+
theme(strip.text.x = element_text(size = 12, face="bold"))+
theme(axis.title.x=element_text(vjust=-0.5))+
scale_y_continuous(breaks = c(0, 0.10, 0.20, 0.30, 0.40, 0.50))+
coord_flip()


## Figure A.9 - display #
	p.A9


# Save figure A.9
  jpeg("FigureA9.jpg")
  p.A9
  dev.off()


################################ TABLE A.10 #################################################################################################################################################################################################
# Variance-covariance  for the individual-level model
load("Estimates.Table1.raw")
Covariance.district.IndividualModel<-rbind(c(apply(Estimates.Table1$Sigma.District,2,mean)[4], HPDinterval(as.mcmc(Estimates.Table1$Sigma.District), p=0.9)[4,]),
c(apply(Estimates.Table1$Sigma.District,2,mean)[1], HPDinterval(as.mcmc(Estimates.Table1$Sigma.District), p=0.9)[1,]),
c(apply(Estimates.Table1$Sigma.District,2,mean)[2], HPDinterval(as.mcmc(Estimates.Table1$Sigma.District), p=0.9)[2,]),
c(mean(Estimates.Table1$Sigma.District[,2]/(sqrt(Estimates.Table1$Sigma.District[,1])*sqrt(Estimates.Table1$Sigma.District[,4])),na.rm=TRUE),
HPDinterval(as.mcmc(Estimates.Table1$Sigma.District[,2]/(sqrt(Estimates.Table1$Sigma.District[,1])*sqrt(Estimates.Table1$Sigma.District[,4]))),p=0.9, na.rm=TRUE)))
Covariance.district.IndividualModel<-rnd(Covariance.district.IndividualModel,2)
rownames(Covariance.district.IndividualModel)<-c("Variance Absenteeism", 
					"Variance - Invalid Voting",
					"Covariance", "Correlation")	
colnames(Covariance.district.IndividualModel)<-c("Mean", "90% CI - Lower", "90% - Upper")


Covariance.election.IndividualModel<-rbind(c(apply(Estimates.Table1$Sigma.Election,2,mean)[4], HPDinterval(as.mcmc(Estimates.Table1$Sigma.Election), p=0.9)[4,]),
c(apply(Estimates.Table1$Sigma.Election,2,mean)[1], HPDinterval(as.mcmc(Estimates.Table1$Sigma.Election), p=0.9)[1,]),
c(apply(Estimates.Table1$Sigma.Election,2,mean)[2], HPDinterval(as.mcmc(Estimates.Table1$Sigma.Election), p=0.9)[2,]),
c(mean(Estimates.Table1$Sigma.Election[,2]/(sqrt(Estimates.Table1$Sigma.Election[,1])*sqrt(Estimates.Table1$Sigma.Election[,4])),na.rm=TRUE),
HPDinterval(as.mcmc(Estimates.Table1$Sigma.Election[,2]/(sqrt(Estimates.Table1$Sigma.Election[,1])*sqrt(Estimates.Table1$Sigma.Election[,4]))),p=0.9, na.rm=TRUE)))
Covariance.election.IndividualModel<-rnd(Covariance.election.IndividualModel,2)
rownames(Covariance.election.IndividualModel)<-c("Variance Absenteeism", 
					"Variance - Invalid Voting",
					"Covariance", "Correlation")	
colnames(Covariance.election.IndividualModel)<-c("Mean", "90% CI - Lower", "90% - Upper")

## # Variance-covariance  for the district-level model
load("Compositional.Cols12.raw")
attach.bugs(Compositional.Cols12, overwrite=TRUE)
Covariance.observation.DistrictModel<-rbind(c(mean(Sigma[,2,2]), HPDinterval(as.mcmc(Sigma[,2,2]), p=0.9)),
c(mean(Sigma[,1,1]), HPDinterval(as.mcmc(Sigma[,1,1]), p=0.9)),
c(mean(Sigma[,1,2]), HPDinterval(as.mcmc(Sigma[,1,2]), p=0.9)),
c(mean(Sigma[,1,2]/(sqrt(Sigma[,1,1])*sqrt(Sigma[,2,2]))),
HPDinterval(as.mcmc(Sigma[,1,2]/(sqrt(Sigma[,1,1])*sqrt(Sigma[,2,2]))))))
Covariance.observation.DistrictModel<-rnd(Covariance.observation.DistrictModel,2)
rownames(Covariance.observation.DistrictModel)<-c("Variance Absenteeism", 
					"Variance - Invalid Voting",
					"Covariance", "Correlation")	
colnames(Covariance.observation.DistrictModel)<-c("Mean", "90% CI - Lower", "90% - Upper")


Covariance.district.DistrictModel<-rbind(c(mean(Sigma.district[,2,2]), HPDinterval(as.mcmc(Sigma.district[,2,2]), p=0.9)),
c(mean(Sigma.district[,1,1]), HPDinterval(as.mcmc(Sigma.district[,1,1]), p=0.9)),
c(mean(Sigma.district[,1,2]), HPDinterval(as.mcmc(Sigma.district[,1,2]), p=0.9)),
c(mean(Sigma.district[,1,2]/(sqrt(Sigma.district[,1,1])*sqrt(Sigma.district[,2,2]))),
HPDinterval(as.mcmc(Sigma.district[,1,2]/(sqrt(Sigma.district[,1,1])*sqrt(Sigma.district[,2,2]))))))
Covariance.district.DistrictModel<-rnd(Covariance.district.DistrictModel,2)
rownames(Covariance.district.DistrictModel)<-c("Variance Absenteeism", 
					"Variance - Invalid Voting",
					"Covariance", "Correlation")	
colnames(Covariance.district.DistrictModel)<-c("Mean", "90% CI - Lower", "90% - Upper")

Covariance.election.DistrictModel<-rbind(c(mean(Sigma.election[,2,2]), HPDinterval(as.mcmc(Sigma.election[,2,2]), p=0.9)),
c(mean(Sigma.election[,1,1]), HPDinterval(as.mcmc(Sigma.election[,1,1]), p=0.9)),
c(mean(Sigma.election[,1,2]), HPDinterval(as.mcmc(Sigma.election[,1,2]), p=0.9)),
c(mean(Sigma.election[,1,2]/(sqrt(Sigma.election[,1,1])*sqrt(Sigma.election[,2,2]))),
HPDinterval(as.mcmc(Sigma.election[,1,2]/(sqrt(Sigma.election[,1,1])*sqrt(Sigma.election[,2,2]))))))
Covariance.election.DistrictModel<-rnd(Covariance.election.DistrictModel,2)
rownames(Covariance.election.DistrictModel)<-c("Variance Absenteeism", 
					"Variance - Invalid Voting",
					"Covariance", "Correlation")	
colnames(Covariance.election.DistrictModel)<-c("Mean", "90% CI - Lower", "90% - Upper")

########################### ESTIMATES TABLE A.10 ###########################
												  #
# Upper panel, Column 2: Observation-specific covariate matrix, 		  #
#	 			  district-level model					  #
Covariance.observation.DistrictModel						  #
												  #
												  #
# Middle panel, Column 1: District-specific covariate matrix, 		  #
#	 			  individual-level model				  #
Covariance.district.IndividualModel							  #
												  #
												  #
# Middle panel, Column 2: District-specific covariate matrix, 		  #
#	 			  district-level model					  #
Covariance.district.DistrictModel							  #
												  #
												  #
# Lower panel, Column 1: Election-specific covariate matrix, 		  #
#	 			 individual-level model					  #
Covariance.election.IndividualModel							  #
												  #
												  #
# Lower panel, Column 2: Election-specific covariate matrix, 		  #
#	 			 district-level model					  #
Covariance.election.DistrictModel							  #
    												  #
###########################################################################

#Backup table in .csv
write.csv(Covariance.observation.DistrictModel, file="TableA10.UpperPanel.csv")
write.csv(Covariance.district.IndividualModel, file="TableA10.MiddlePanelLeft.csv")
write.csv(Covariance.district.DistrictModel, file="TableA10.MiddlePanelRight.csv")
write.csv(Covariance.election.IndividualModel, file="TableA10.LowerPanelLeft.csv")
write.csv(Covariance.election.DistrictModel, file="TableA10.LowerPanelRight.csv")



